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Chapter 2 


Time 


Time is the primary independent variable in GMAT. 
Time is used in integrating the equations of motion, and 
calculating planetary ephemerides, the orientations of plan- 
ets and moons, and atmospheric density among others. 
GMAT uses three types of time systems depending on 
the type of calculations being performed: universal time 
systems based on the Earth’s rotation with respect to the 
Sun: dynamic time systems that are based on the dy- 
namic motion of the solar system and take into account 
relativistic effects; and atomic time systems based on the 
oscillation of the cesium atom. Each of these time sys- 
tems has specific uses and is discussed below. In addition, 
universal, dynamic, and atomic time systems can be ex- 
pressed in different time formats. The two time formats 
used in GMAT are the Modified Julian Date (MJD) for- 
mat, and the Gregorian Date (GD) format. In the next 
section, we’ll take a look at the time systems used in 
GMAT, and when GMAT uses each time system. Then 
we’ll look at the different time formats. 


2.1 Time Systems 


GMAT uses several different time systems in physical 
models and spacecraft dynamics modelling. The choice 
of time system for a particular calculation is determined 
by which time system is most natural and convenient, 
as well as the accuracy required. In general, for deter- 
mining Earth’s orientation at a given epoch, we use one 
of several forms of Universal 'lime (UT), because uni- 
versal time is based on the Earth’s rotation with respect 
to the Sun. Planetary ephemerides are usually provided 
with time in a dynamic time system, because dynamic 
time is the independent variable in the dynamic theories 
and ephemerides. The independent variable in spacecraft 
equations of motion in GMAT is time expressed in an 
atomic time system. Let’s look at each of these three 
systems, starting with atomic time. 


2.1.1 Atomic Time: TAI and A.l 

Atomic time (AT) is a. highly accurate time system which 
is independent of the rotation of the Earth. 1 Therefore, 
AT is a natural system for integrating a spacecraft’s equa- 
tions of motion. AT is defined in terms of the oscillations 
of the cesium atom at mean sea, level. The duration of 
the SI second is defined to be 9,192,631,770 oscillations of 
the cesium nuclide 133Ce. Two atomic times systems are 
are used in GMAT: A.l, and international atomic time 
(TAI). A.l is in advance of TAI by 0.0343817 seconds. 

A.l = TAI + 0.0343817sec (2.1) 

where 

TAI — UTG + A AT (2.2) 

and A/17’ is the number of leap seconds, added since 1972, 
needed to keep i UTG— UT1 |< O.dsec. GMAT reads 
A AT from the file named tai-v.tc.dab. For times that ap- 
pear before the first epoch on the file, GMAT uses the 
first value found in the file. For times that appear after 
the last epoch, GMAT uses the last value contained in the 
file. Currently, GMAT uses A.l time as the independent 
variable in the equations of motion. TAI is used as a time 
system for defining spacecraft state information. 

Now let’s look at the universal time system. 

2.1.2 Universal Time: UTC and UT1 

All of the universal time (UT) scales are based on the 
Earth’s rotation with respect to a fixed point (sidereal 
time) or with respect to the Sun (solar time). The ob- 
served universal time (UTO) is determined from observa- 
tions of stellar transits to determine mean local sidereal 
time. UT1 is UTO corrected for the Earth’s polar mo- 
tion and is used when the instantaneous orientation of the 
Earth is needed. UTC is the basis for all civil time stan- 
dards. It is also known as Greenwich mean time (GMT) 
and Zulu time (Z). The UTC time unit is defined to be 
an SI second, hut UTC is kept within 0.9 seconds of UT1 
by occasional leap second adjustments. The equation re- 
lating UTC and UT1 is 

UTC — UT1 — AUT1 (2.3) 
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In GMAT, AUT1 is read from the file eopcO/. 62-now pro- 
vided by the International Earth Rotation and Reference 
Systems Service (IERS). The file containing the latest 
measurements and predictions can be found at 
http://www.iers.org/. For times past the last epoch con- 
tained in the file. GMAT uses the last value of Af/Tl 
contained in the file. GMAT uses UTC as a time system 
to define spacecraft state information. UT1 is used to 
determine the Greenwich hour angle and for the sidereal 
time portion of FK5 reduction. 


2.1.3 Dynamic Time: TT, TDB and TCB 

Dynamical time is the independent variable in the dy- 
namical theories and ephemerides. This class of time 
scales contains terrestrial time (TT), Barycentric Dynam- 
ical time (TDB), and Barycentric Coordinate Time (TCB). 
TDB is the independent variable in the equations of mo- 
tion referred to the solar system barycenter. It is also the 
coordinate time in the theory of general relativity. De- 
spite the fact that the Jet Propulsion Laboratory (JPL) 
J2000.0 ephemerides are referred to in TDB, TT is fre- 
quently used. This is because TT and TDB always differ 
by less than 0.002 sec. As higher accuracy or more sensi- 
tive missions are planned, the difference may need to be 
distinguished. In this section we’ll discuss how to calcu- 
late TT, TDB and TCB, and discuss where each is used 
in GMAT. 

TT is the independent variable in the equations of 
motion referred to the Earth’s center. It is also the proper 
time in the theory of general relativity. The unit of TT is 
a day of 86400 SI seconds at mean sea level. In GMAT, 
TT is used in FK5 reduction, and as an intermediate time 
system in the calculation of TDB and TCB. TT can be 
calculated from the following equation: 

TT = TAl + 32.184 sec (2.4) 

Calculating TDB exactly is a complicated process that 
involves iteratively solving a transcendental equation. For 
this reason, it is convenient to use the following approxi- 
mation. 

TDB « TT + 0-001658 sin + 0.00001 385 sin 2M E 

V ^ -- 

units of seconds 

(2.5) 

Note that the term in the underbrace has units of seconds, 
and depending upon the units of TT, which is usually 
in days, a conversion of the term may be necessary be- 
fore performing the addition with TT. Me is the Earth’s 
mean anomaly with respect to the sun and is given ap- 
proximately as 

M e « 357.5277233 + 35, 999.050347 tt (2.6) 


where 'T'tt is the time in TT expressed in the Julian Cen- 
tury format. T tt can be calculated from 

Ttt = J J>™ Z*>**h™* 

11 ■-|,j rqr \ ' / 

-> 0 , o/o 

where J D-tt is the time in TT expressed in the Julian 
Date format. For a more complete discussion of the TDB 
time system, see Vallado 1 (pp. 195-198) and Seidelmann 2 
(pp. 41-48). GMAT uses TDB as the default time system 
in the JPL ephemerides files. There is an option to use TT 
in the ephemerides Using the UseTTForEphemeris flag. 

The last dynamic time system GMAT uses is Barycen- 
tic Coordinate Time. In 1 992, the IATJ adopted this sys- 
tem and clarified the relationships between space-time co- 
ordinates. 2 In general, calculating TCB requires a four- 
dimensional space-time transformation that is well be- 
yond the scope of this discussion. However, TCB can be 
approximated using the following equation: 

TCB - TDB = L b (JD - 2443144.5)86400 (2.8) 

The present estimate of the value of Lb is 1.550505x10 “ 8 
(+/ — 1x10“ ,• (Fukushima et al., Celestial Mechanics, 
38, 215, 1986). Tt is important to note that the main dif- 
ference between TDB and TCB is a, secular drift, and that 
as of the J2000 Epoch, the difference was approximately 
11.25 seconds and grooving. 3 GMAT uses time in the 
TCB system to evaluate the IAU data for the spin axes 
and prime meridian locations of all planets and moons ex- 
cept for Earth. 4 Note that Seidelmann 4 mistakenly says 
that time in TCB should be used in the equations given 
for the pole and meridian locations of the planets. The 
correct time to rase is TDB, and GMAT uses this time 
system . 


2.2 Time Formats 

There are two time formats that GMAT uses to repre- 
sent time in the systems discussed above. These formats 
are called the Gregorian Date (GD), and the Julian Date 
(JD). The difference between the GD and JD formats is 
how they represent the Year, Month, Day, Hours, Min- 
utes, and Seconds of a given date. The GD format is 
well known, and the J2000 epoch is expressed as, 01 Jan 
2000 12:00:00.000 TT. The reference epoch for the GD 
calendar is the beginning of the Christian Epoch. The 
JD format represents an epoch as a continuous number 
containing the day and the fraction of day. 

The J2000 epoch is commonly used in astrodynamics 
as a reference epoch for planetary and other data. The 
J2000 epoch occurred at 01 Jan 2000 12:00:00.000 TT. 
The time system, TT, is important for precise applica- 
tions! While the J2000 epoch is a specific instant in time, 
the numerical value changes depending upon which time 
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system you express it in. We can make an analogy with 
vector algebra where we have an abstract quantity that 
is a vector, and can’t write down a set of numbers rep- 
resenting the vector until we choose a coordinate system. 
Similarly, the J2000 epoch can be written in any of the 
different time systems and formats. All of the following 
are equivalent definitions of tire J2000 Epoch: 


minutes, and seconds respectively. This equation is valid 
for the time period 01 Mar. 1900 to 28 Feb. 2100. 

For numerical reasons it is often convenient to work in 
a Modified Julian Date (MJD) format to ensure we can 
capture enough significant figures using double precision 
computers. In GMAT the MJD system is defined as 


2451545.0 TT 
2451544.9996274998411 TAI 
2451 544. 9992571 294708 UTC 
2451544.9999999990686 TUB 

01 Jan 2000 12 : 00 : 00.000 TT 
01 Jan 2000 11 : 59 : 27.815986276 TAI 
01 Jan 2000 11:59:55.815986276 UTC 
01 Jan 2000 11 : 59 : 59.999919534 TDB 


In the next two sections we’ll look at how to convert 
an epoch in a given time system from the GD format to 
the JD format, and vice versa. 


MJD = J D — 2, 430, 000.0 (2.10) 

where the reference epoch expressed in the GD format is 
05 Jan 1941 12:00:00.000. However, -we must be careful 
in calculating the Modified Julian Date, or we will lose 
the precision we are trying to gain. GMAT calculates the 
MJD as follows: 


JDay =367 K — lnt 


( 7(Y + lnt 


'M+9\\\ 


12 


V 


+ 


J 


( 27oM \ 

hit ( - j + D + 1,721, 013.5 - 2, 430, 000.0 

( 2 . 11 ) 


2.2.1 Julian Date and Modified Julian Date 


PartofDay = 


s 

( -m 

60 

60 


+ U 


24 


( 2 . 12 ) 


The Julian date is a time format in which we can express a 
time known in any of the Atomic, Universal or Dynamic 
time systems. The Julian Date is composed of the Ju- 
lian day number and the decimal fraction of the current 
day. SeidelmamT (pp. 55-56) says “The Julian day num- 
ber represents the number of days that has elapsed, at 
Greenwich noon on the day designated, since ...the epoch 
noon Jan 1 4713 B.C. in the Julian proleptic calendar. 
The Julian date (JD) corresponding to any instant is, by 
simple extension to this concept, the Julian day number 
followed by the fraction of the day elapsed since the pre- 
ceding noon” . 

The fundamental epoch for most astrodynamie calcu- 
lations is the J2000.0 epoch. 2 This epoch is GD 01 Jan 
2000 12:00:00.000 in the TT time system and is expressed 
as JD 2451545.0 TT. To convert between Julian Date for- 
mat and Gregorian Date format, GMAT uses Algorithm 
14 from Vallado 1 


MJD = (JDay - 2, 430, 000.0) + PartofDay (2,13) 

The important subtlety is that we must subtract the MJD 
reference from the JD, before we add the fraction of day, 
to avoid losing' precision in the MJD. 


2.2.2 Gregorian Date 

The Gregorian Date format is primarily used as a time 
system in which to enter state information in GMAT. GD 
is not a convenient time format for most mathematical 
calculations. Hence, GMAT often takes input in the GD 
format and converts it to a MJD format for use internally. 
The algorithm for converting from GD to MJD is taken 
from Vallado 1 and reproduced here verbatim. 


A 


JD = 367 Y - lnt 


7 Y + hit 


M + 9\\ \ 


12 


lnt 


m 


60 


+ m 


+ D+ 1,721,013.5 + 


60 


f 

/ (2.9) 

+ H 


24 


where Y is the four digit year, lnt signifies real trun- 
cation, and Ad , D, H, m and s are month, day, hour, 


JD — 2, 415,019.5 
; 190(1 ~ 36A25 

Year = 1900 + TRUNC(T; 9 (m) 

LeapYrs = TRUNC((K ear — 1900 — 1 j (0.25)) 

Days = (JD - 2,415, 019.5)- 

(( Year - 1900)(365.0) + LeapYrs ) 

IF Days < 1.0 THEN 



14 


CHAPTER 2. TIME 


Y ear — Y ear — 1 

LeapYrs — TRTJNC((V' eai — 1900 — 1)(0.25)) 

Days — (JD 2, 415, 019.5)— 

((Year — 1900) (365.0) + LeapYrs) 

If ( Year Mod 4) = 0 Then 

LMonth[ 2] — 29 
DayofYr = TRUNC (Days) 


Sum days in each month until 

LMonth + 1 summation > DayofYr 

M on, — of months in summation 

Day = Dayof Yr — LM orttli summation 

r — ( Days — DayofYr) 24 

h = TRUNC(Temp) 


min - T RUN C ( ( 7 ’ ernp - h) 60) 

a = ('Temp - h. - ~ N ) 3600 
V. 60 j 

2.3 Conclusions 


In this chapter we looked at three time systems that 
GMAT uses to perform internal calculations: atomic time, 
universal time, and dynamic time. Atomic time is used 
to integrate spacecraft equations of motion, while univer- 
sal time is used to determine the sidereal time and green- 
which hour angle for use in FK5 reduction. Dynamic time 
systems are used in the JPL ephemerides and in the IAU 
planetary orientation data. Time, in any of these time 
systems, axe represented in two formats: the Gregorian 
Date, and Julian Date. We looked at how to convert be- 
tween different time systems, and between different time 
formats. 



Chapter 3 

Coordinate Systems 


There are numerous coordinate systems used in space 
mission analysis, that when used appropriately can greatly 
simplify the work and yield insight that is not obvious 
otherwise. Some examples are equatorial and ecliptic sys- 
tems, and rotating coordinate systems based on the rel- 
ative motion of two bodies such as the Earth and Moon. 
GMAT has the capability to calculate many parameters 
in different coordinate systems, and these parameters can 
then be used in plots, reports, solvers, control flow state- 
ments and stopping conditions to name a few. 

In this chapter we investigate how GMAT performs 
coordinate system transformations, and how different co- 
ordinate systems are defined. We begin by defining some 
notation. Next, we look at how' to transform a vector and 
its first derivative from one coordinate system to another 
when the coordinate systems are translating and rotating 
with respect to one another. Finally, we look at each co- 
ordinate system defined in GMAT and discuss how to find 
its rotation matrix and the first derivative of the rotation 
matrix to rotate to the J2000 coordinate system. 


3.1 General Coordinate System 
Transformations 


GMAT has the capability to take a position and veloc- 
ity vector in one coordinate system, and convert them to 
another coordinate system that may be both translating 
and rotating with respect to the original system. In this 
section we derive the equations governing coordinate sys- 
tem transformations and describe the algorithm GMAT 
uses to transform position and velocity vectors. 

Let’s start by defining some notation. In Fig. 3.1, 
we see an illustration of a point, o drawn with respect to 
two coordinate systems and Tf. The vector r° is the 
position vector of point o with respect to frame Ti. The 
vector Tj is the position vector of point o with respect to 
Tf. r fi is the vector from the origin of Ti to T < . Let’s 
define the rotation matrix that rotates from Ti to Tf as 
R fi. Finally, let’s define the angular velocity u>f,, as the 
angular velocity of Ti with respect to Tf. To simplify 


the notation, we assume that a vector is expressed in the 
frame denoted by the right-most subscript. If we need 
to express a vector in another coordinate system, we use 
curly brackets. As an example, {r?}/ is the position of 
point o, with respect to frame Ti, expressed in frame Tf. 
In summary, we have 


r? 

Position vector of point o w/r/t frame T; 
expressed in T, 

K}/ 

Position vector of point o w/r/t frame T, 
expressed in Tf 

r/,: 

Position vector from origin of Ti to origin 
of Tf expressed in Tf 

R fi 

Rotation matrix from frame Ti to Tf 

Ufi 

Angular velocity of frame Ti w/r/t Tf, 
expressed in frame Ti 

Wfi}f 

Angular velocity of frame T w/r/t Tf, 
expressed in frame T / 


To further simplify the notation, let’s drop the su- 
perscript “o” from r? and rj. Then, from inspection of 
Figure 3.1 we can write 


Tf = n.f i T i + tif 
Rot. Trails. 


(3.1) 


Equation (3.1) is the equation used to convert a vector 
known in frame Ti to a vector in frame Tf , where both a 
rotation and a translation are required. The first term in 
Eq, (3.1) is the term that performs the rotation portion of 
the transformation. Here, n is the position vector w/r/t 
to Ti and is expressed in T. R/, is the rotation matrix 
that rotates from T.- t to Tf. t if is the vector that goes 
from the origin of Tf to the origin of T ,, and is expressed 
in Tf. 

We also need to be able to determine the time rate 
of change of a vector in frame Tf if we know the time 
rate of change of the vector in T; t . To determine the 
equation that describes the transformation, we must take 
the derivative of Eq. (3.1) with respect to time. 

dr f dRfiTi , dr if 

dt ~ di " + df 1 


15 
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Figure 3.1: Illustration of a Translating and Rotating Coordinate System 


Let’s use a single dot above a variable to denote the first 
derivative of that variable with respect to time. Then, we 
can expand this to obtain 

if — R fiTi + Rfiti + r,-/ (3.3) 

In Eq. (3.3) we see a term that contains the time deriva- 
tive of the rotation matrix from T to Tf. We can write 
the time derivative of as 

R/i = R/i w /a = { w /i)/R/i (3-4) 

where w/* is the angular velocity of Ti with respect to 
Tf expressed in T. The skew symmetric matrix, u> x , is 
defined as 

( 0 ~‘"'Z Uy \ 

c j z 0 — J (3.5) 

' 0 J 

In summary, using Eq. (3.4) to transform a derivative 
vector from Ti to Tf we can use any of the following 
three equations: 

if = R/,uiy,r, + Rfiii + r</ (3.6) 

Hot. Trans. 

if = {titfijfRfiTi + Rfiii + iif (3.7) 

Frans. 


if - RfiTj + Rfiij + iif^ (3.8) 

Rot. T rans. 

We choose between Eqs. (3.6), (3.7), or (3.8) depending 
on the type of information we have, and which frame is 
most convenient, to express the angular velocity u> fi in. 
In general, we know r< and q. To perform the transfor- 
mation we need to determine R, R, and iif and these 
quantities depend on Ti and Tf. 

One of the difficulties in implementing coordinate sys- 
tem transformations in GMAT is that we often can’t cal- 
culate R fi and Rj, directly. For example, it is nontrivial 
to directly calculate the rotation matrix from the Earth 
fixed frame to the Moon fixed frame. Hence, we need, to 
choose a convenient intermediate coordinate system. We 
choose the axis system defined by Earth’s mean equinox 
and mean equator at the J2000 epoch, denoted Tj 2k , as 
the intermediate reference frame for all transformations 
that require an intermediate transformation. This choice 
is motivated by the fact that most of the data needed to 
calculate R and R is given so that it is fast and convenient 
to calculate R j 2 k ,i, and R .j ik .i. 

The steps taken to perform a general coordinate trans- 
formation in GMAT are described below and illustrated 
in Fig. 3.2. We start with a vector and its first derivative 
known in frame T, and wish to determine the vector and 
its first derivative with respect to frame Tf. However, we 


Rot. 
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assume that the transformation to go directly from Ti to 
Tf is not known. 

The first step in the process is to perform a rotation 
from Ti to Tj. ih ■ We define this intermediate system as 
Tf. No translation is performed in step one. Using only 
the rotation portions of from Eqs. (3.1) and (3.8) we see 
that 

(r<}i = Rjjfc.iH (3.9) 

{ U } | ll./ 2i .. .jr, T (3. iO) 

The second step is to perform a translation from the origin 
of Ti to the origin of Tf. We define this second interme- 
diate system as Tv- Tv has the same origin as Tf but has 
the same axes as Tj. ik . From inspection of Fig. 3. 2 we can 
see that 

foil = + {rfii}j 2k + {r /} 2 (3.11) 

Solving for r j we obtain 

{t;} 2 - {rOx - {rRi). hk - (3-12) 

where {r Ri } j is the vector from the origin of Ti to the 
origin of Tr expressed in Ti 2 k- Similarly {r }R.}j 2k is the 
vector from the origin of Tr to the origin of Tf expressed 
in Tjvk • Because the vector {r /} 2 is expressed in an iner- 
tial system we can we can take the derivative of Eq .(3.12) 
to obtain 

{ V/ } 2 - {v^ - {vRi}. hk - {v/fi) j2 . (3.13) 

where {vr<} j is the velocity of the origin of Tr w/r/t 
the origin of Ti. Similarly. { v fR.}j 2k is the velocity of the 
origin of Tf w/r/t the origin of Tr. Finally, we perform a 
rotation from Tj 2k to Tf about the origin of Tf to obtain 
the desired quantities. 

r f = (v/} 2 (3.14) 

if = R f .. hk {r /} 2 + R f ,j M {v/} 2 (3.15) 


3.2 Pseudo-Rotating Coordinate 
Systems 


In mission analysis, sometimes it is useful to consider a ro- 
tating coordinate system to be inertial at. a given instant 
in time. In this case, we ignore the effects of rotation 
on the velocity. Let’s call systems where we neglect the 
rotational effects on velocity pseudo-rotating coordinate 
systems. 

To perform transformations to a pseudo-rotating coor- 
dinate system, the equations to convert, a position vector 
do not change and art: given by Eqs. (3.9) and (3.14). 


However, the velocity conversion equations change be- 
cause we neglect the terms that contain R. For pseudo- 
rotating coordinate systems the velocity transformations 
shown in Eqs. (3.10) and (3.15) become 


dr. : 

dt 


— Rjafc.i 

1 


dn 

at 


(3.16) 


and 


dvf 

dt 


= R S.Hk { v /} 2 


(3.17) 


To perform the transformations describe in the last 
few sections, we need to be able to calculate the rotation 
matrix between any coordinate system and Tj ?k , and the 
derivative of the rotation matrix. In the following sec- 
tions we calculate these matrices for the systems used 
in GMAT. We assume that we want the transformation 
from some generic frame T; to Tj. lk which is the Earth 
Mean J2000 Equatorial (MJ2000Eq) system. The rota- 
tion matrix from Tj 2k to Ti can be found by the simple 
relationship. 


and 


R-vh* = R/i-i = RL,-: (3.18) 

Rtf = R jLi = (3-19) 


3.3 The Tj Vk Inertial System and 
FK5 Reduction 


It is well know that Newton’s laws must be applied in an 
inertial system. The struggle to determine a truly inertial 
system has continued since Newton’s time. In reality, the 
best we can do is approximate a truly inertial system in 
which to apply Newton’s laws. In GMAT that system is 
the FK5 system, here called Tj ik . The Tj 2k system is 
referenced to the Earth’s equator and the Earth’s orbit 
about the sun. Because neither of these two planes are 
fixed in space, we must pick an epoch and define an in- 
ertial system based on the geometry at that epoch. This 
epoch is commonly chosen as the J2000 epoch. In this 
section, we present the definition of the Tj 2k system, and 
discuss the transformation from Tj 2k to the Earth Fixed 
system. This transformation is called FK5 reduction. We 
begin with a conceptual discussion of how the Earth’s 
spin axis moves with respect to inertial space. We con- 
clude this section with a presentation of the mathematical 
theory of FK5 reduction. 


3.3.1 Overview of FK5 Reduction 


The inertial system most commonly used in astrodynam- 
ics as of this writing is the FK5 system. We call this 
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system Fjn- The .Tj 2k system is used for many calcu- 
lations in GMAT. The two most important are for inte- 
grating equations of motion, and as an intermediate sys- 
tem for coordinate system transformation. Tj„ k is used 
throughout the astrodynamics community as the coordi- 
nate system to represent time varying data such as plane- 
tary ephexnerides and planetary pole and prime meridian 
locations. 

The rigorous mathematical definition of Tj. zk is com- 
plex. So, let’s start with a simple qualitative explana- 
tion. The nominal 2 -axis of Tj 2k is normal to the Earth’s 
equatorial plane. The nominal x-axis points along the 
line formed by the intersection of the Earth’s equatorial 
plane and the ecliptic plane, in the direction of Aries. 
The nominal y-axis completes the right-handed system. 
Both the equatorial and ecliptic planes move slowly with 
respect to inertial space. The rigorous definition of FK5 
uses the mean planes of the ecliptic and equator, at the 
J2000 epoch. We’ll take a closer look at the mathematical 
definitions of the mean ecliptic and equator below. 

FK5 reduction is the transformation that rotates a 
vector expressed in the T system to the Earth Fixed 
coordinate system. To perform this transformation obvi- 
ously requires an understanding of how the Earth’s ori- 
entation changes with respect to Tj . ik . The time varying 
orientation of the Earth is complex and is due to com- 
plicated interactions between the Earth and it's external 
environment and complicated internal dynamics. In fact, 
the dynamic orientation of the Earth is so complicated 
that we can’t model it completely and FK5 reduction is 


a combination of dynamics models and empirical obser- 
vations that are updated daily. 

We describe the orientation of the Earth using three 
types of motion. The first type, including precession and 
nutation, describes how the Earth’s principal moment of 
inertia moves with respect to inertial space. 2 The motion 
is illustrated in Fig. 3.3. The principal moment of inertia 
is defined as the Celestial Ephemeris Pole, and due to the 
fact that Earth’s mass distribution changes with time, the 
Celestial Ephemeris Pole is not constant with respect to 
the Earth’s surface. Precession is the coning motion that 
the Celestial Ephemeris Pole makes around the ecliptic 
north pole. The other principal component, of the mo- 
tion of the Celestial Ephemeris Pole is commonly called 
nutation and is the oscillation in the angle between the 
Celestial Ephemeris Pole and the north ecliptic pole. The 
theory of Precession and Nutation come from dynamics 
models of the Earth’s motion. The second type of motion 
is called sidereal time, and represents a rotation about the 
Celestial Ephemeris Pole. The sidereal time model is a 
combination of theory ad observation. The third motion 
is that of the Earth’s instantaneous spin axis with respect 
to the Earth’s surface. As, we’ll see below, the Earth’s 
spin axis is not. constant with respect to the Earth’s crust 
and it’s motion is called Polar Motion. A portion of polar 
motion is due to complicated dynamics, and a portion is 
due to unmodelled errors in nutation. Polar motion is 
determined from observation. Now that we’ve had a brief 
introduction to precession, nutation, sidereal time, and 
polar motion, let’s look at each in more detail. 



Ecliptic North Pole 



Figure 3.3: Inertial Motion of Earth’s Spin Axis 


Sidereal Time 

3.3.2 Precession Calculations 



J L)tdb & JDtt 

(3.20) 


JV-tdb- 2,451, 

545.0 . 

1'tdb 

36525 

13.21) 

2306.2181” Tt-db + 0.301882£ DB 

+ 0.0179982 % DB 



(3.22) 

2004.3109 "T t 

db ~ 0.42665'i? DB 

- 0.041 8332 % DB 



(3.23) 

2306.2181 "T TD b + 1.094682?. BB 

+ 0.0182032:^,3 



(3.24) 

P = 

R,(- 2 )R 2 (©)R 3 (- 

-0 (3.25) 


3.3.3 Nutation Calculations 

GMAT has the ability to use either the 1980 IAU Theory 
of Nutation, or the IERS 1996 Theory of Nutation. There 
are some calculations that are common to both, so let’s 
look at. them first. The mean obliquity of the ecliptic at 
the J2000 epoch, e, is given by 
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Figure 3.4: Intermediate Transformations and Coordinate Systems in FK5 Reduction 


e —84381.448 - 46.81507r DB 


+ 0.0018137 


■3 

TDB 


0.00059 1% B 


(3.26) 


n — 1 25.04455501° + (-6962890.26657 tdb+ 

7.47227 £ db + 0.0077027 £ DB - 0.000059397^)" 

(3.31) 


As we mentioned previously, Earth’s nutation is caused 
by the combined gravitational effect of the Moon and Sun. 
So, we would expect to see the time dependent location of 
the Moon and Sun appear in the equations for Earth nu- 
tation. The theories of nutation described below take into 
account of the Moon and Sun position by modelling mean 
anomalies of the Moon and Sun, l and V respectively, the 
mean argument of latitude of the Moon, F, the difference 
between the mean longitude of the Sun and Moon, D, and 
the mean longitude of the ascending node of the Moon’s 
orbit, f 1. The equations used to determine these values 
as a function of 1'tdb are: 

l —134.96340251° + (1717915923. 21787tdb+ 

31. 87927 Ip DB + 0.051 6357 £ db - 0.000244707 £ BB )" 

(3.27) 

l ' =357.52910918°+ (129596581.04817tu.u- 

0. 55327 £ bb - 0.0001367 y DB - 0.000011497^)" 

(3.28) 

F =93.27209062° + (1739527262.84787tx.b- 

12.75127$ bb + 0.001 0377 y BB + 0.Q00004177~1 DB )” 

(3.29) 

D =297.85019547°+ (1602961601. 20907 tur- 

6.37067f OB + 0.0065937’| £)B - 0.000031697^^” 

(3.30) 


1980 Nutation Theory 


106 

AS&iqso — (Ai + BiTtdb ) sin Up 

1=1 

(3.32) 

100 

Acisso = ^2 ( Ci + 71/7 tdb) cos dp 

i - 1 

(3.33) 

a v ----- aul + anV + a^F + 0,4 iD + 05 * fi 

(3.34) 

N = R 1 (-e)R 3 (-AT)R 1 (e) 

(3.35) 


1996 Nutation Theory 

The 1996 theory of nutation published by the IERS is a 
higher fidelity model of Earth nutation. There are two 
primary differences between the 1908 IAU theory and the 
1996 IERS theory. The first difference is the 1996 the- 
ory uses a 263 term series expansion for the effects of the 
Earth and Moon. The 1980 theory uses a 106 term series. 
The second difference is that the 1996 theory has a sec- 
ond series expansion to account for the effects of nutation 
caused by the more massive planets. The planetary series 
expansion is a 118 term series. Let’s begin with the equa- 
tions for the Earth and Moon’s effect on Earth nutation, 
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according to the 1996 IERS theory: 

263 

A^1996 — 52 + BiT'rDB ) sin Op + Ei cosa P (3.36) 

i—1 

263 

A6i996 = 5_, + E/Ttob) cosa p + Ei sina p (3.37) 

i = 1 

<i p — auMc + avMo + asfjUjtfo + a^/J© + (3.38) 

To calculate the planetary effects on nutation, we be- 
gin by calculating the mean heliocentric longitude of the 
planets. Only the effects of Venus(V), Mars(M), Jupiter J), 
and Saturn(S) are included in the theory. We require the 
Earth’s (E) mean longitude also. The mean longitudes 
are calculated using: 

A v = 181.979800853° + 58, 517.8156748 r ir£>B 
A E - 100.466448494° + 35, 999.3728521'/t£>b 
X M = 355.433274605° + 19, 140.299314 J' TDB 
A j = 34.351483900° + 3, 034.905674647 tdb 
A s = 50.0774713998° + l,222.11379404'i’ T£ , B 

The general precession in longitude, p a , is calculated us- 
ing 

p a - 1 .396971 37214°7 tdb + 0.00030867 tdb 


3.3.4 Sidereal Time Calculations 

To calculate the sidereal time of the Earth, we need the 
current time which is then used to determine the Green- 
wich Mean Sidereal Time (GMST) and the equation of 
the equinoxes. GMST is calculated using: 

Ogmst =1 .0096582261 5e6 + 4. 74660027721 9299el07;/ ri 
+ 1. 3965607 $ r , + 9.3e — 57f, ri (arcseconds) 

(3.46) 

The calculation of the equation of the equinoxes is de- 
pendent upon the time. If the Julian date falls after 
2450449.5, then we use 

EQequinox = A 4' cos £ + 0.00264 sin Q + 0.000063’ sin 211 

(3.47) 

If the Julian date falls on or before 2450449.5 we use 


EQequinox ~ z: COS C 

(3.48) 

0 A$T = @GM ST + EQequinox 

(3.49) 

ST = R'3 

(3.50) 

Polar Motion Calculations 


>, 3 

§ 

1! 

P3 

to 

T 

P 

i— * 

T 

5^ 

(3.51) 


Finally, the planetary terms are calculated using: 

118 3.4 Deriving Rj 2fcji and R j 2k>i for 

— 52 (At + BiTtdb) sincipi (3.39) Various Coordinate Systems 

i—1 


118 

A e p i — 52 (C + D/1tdb) cos a p i (3.40) 

i=l 

a p i = anXv + «2Ab + «3»A m + a.uXj + a-a X$ + 

«6 iVa + 0,-jiD + a^F + a.Qil + am.jQ (3.41) 

Aty = A^iggfl + A \|V (3-42) 

optional 

Ae = AeiQ96 + A e p i (3.43) 

optional 

e = e + Ae (3.44) 

In GMAT, the planetary terms are optional. If the user 
has selected to include the planetary terms, the 


In GMAT, there are numerous coordinate systems that 
can be used to define variables, s topping conditions, space- 
craft. states, and other quantities. Some examples include 
the Earth centered mean ecliptic of .12000 system, the 
Earth-fixed system, the Mars equator system, and the 
Earth-Moon rotating system. 

In the following subsections, we determine how R .> 2k . t 
and Rq 2fci j are calculated in GMAT for all of the coordi- 
nate systems available in GMAT. Let’s begin by looking 
at coordinate systems defined by the equator of a celestial 
body. 

3.4.1 Equator System 

The Equator axis system has the following nominal con- 
figuration: 

• z-axis: Along the line formed by the intersection of 
the bodies equator arid the ecliptic plane. 


N = R.i(— e)R 3 (— A40Ri(e) 


(3.45) 
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Figure 3.5: IAU Definition of Pole and Meridian Locations for Planets and Moons 


t„ is a constant value for the TOEEq system. For a given 
t Q , the matrices associated with the TOEEq system only- 
need to be evaluated once and can be reused later when 
necessary. P(f„) and N(f 0 ) are part of the FK5 reduc- 
tion algorithm and are explained in detail in Sec. 3.3.3 
and 3.3.2 . Because t 0 is fixed for a TOEEq system, the 
derivative of R/< is identically equal to the zero matrix. 

/0.0 0.0 o.o\ 

R = 0.0 0.0 0.0 (3.61) 

\ 0.0 0.0 0 . 0 / 

3.4.4 Mean of Epoch Equator (MOEEq) 

The Mean of Epoch Equator axis system is defined as 
follows: 

• cc-axis: Along the mean equinox at the chosen epoch. 
The axis points in the direction of the first point of 
Aries. 

• y-axis: Completes the right-handed set. 

• 2 -axis: Normal to the Earth’s mean equatorial plane 
at the chosen Epoch. 

The MOEEq is an intermediate system in FK5 reduction 
and R K and R fi for the MOEEq system can be calculated 
using the following equations 

R/i = P r (t 0 ) (3.62) 

where t a is the epoch defined in the coordinate system 
description provided by the user in the epoch field. Hence 
t 0 is a constant value for the MOEEq system. For a given 
t 0 , the matrices associated with the MOEEq system only 
need to be evaluated once and can be reused later when 
necessary. P(t 0 ) is described in Sec. 3.3.2. Because t 0 is 


fixed for a MOEEq system, the derivative of R/; is the 
zero matrix. 

/ 0.0 0.0 0 . 0 \ 

R = 0.0 0.0 0.0 (3.63) 

\ 0.0 0.0 0 . 0 / 

3.4.5 True of Date Equator (TODEq) 

R u and R/j for the TODEq system can be calculated 
using the following equations Valiado 1 Fig. 3-29). 

R/i = N r (f e )P T (t,) (3.64) 

where t a is the epoch. Unlike the TOEEq sytem, for the 
TODEq system t a is a variable and usually comes from 
the epoch of the object whose state we wish to convert. 
P (t Q ) and N(to) are? part of the: FK5 reduction algorithm 
and can be found in Valiado pgs. 214 - 219. Because t 0 
is not fixed for a TODEq system the derivative of R H 
is non-zero. However, we will assume its derivative is 
negligibly small so that 

/ 0.0 0.0 0 . 0 \ 

R= 0.0 0.0 0.0 (3.65) 

\ 0.0 0.0 0 . 0 / 

3.4.6 Mean of Date Equator (MODEq) 

R/i and R/i for the? MODEq system can be? calculated 
using the following equations 

R/i - P T (t„) (3.66) 

where t a is the epoch. Unlike the MOEEq sytem, for the? 
MODEq system t 0 is a variable and usually comes from 
the epoch of the object whose state we wish to convert. 
P(t 0 ) and N(f 0 ) are: part, of the? FKo reduction algorithm 
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and can be found in Vallado 1 pgs. 214 - 219. Because 
t 0 is not fixed for a MODEq system, the derivative of 
R i i is non-zero. However, we will assume its derivative 
is negligibly small so that 

/0.0 0.0 o.o\ 

R = 0.0 0.0 0.0 (3.67) 

\ 0.0 0.0 0 . 0 / 

3.4.7 Mean of Date Ecliptic (MODEc) 

Rji and R/j for the MODEc system can be calculated 
using the following equations 

R/i = P r (t 0 )R?'(e) (3.68) 

where t 0 is the epoch. For the MODEc system t a is a 
variable and usually comes from the epoch of the object 
whose state we wish to convert. P (t a ) comes from the 
FK5 reduction algorithm and can be found in Vallado 1 
pgs. 214 - 2.19. e is given by Vallado, 1 Eq. (3-52). For 
a more complete discussion, you can also refer to Seidel- 
mann, 2 pgs. 114 - 115. Because t a is not feted for a 
MODEc system, the derivative of is non-zero. How- 
ever, we will assume its derivative is negligibly small so 
that 

/ 0.0 0.0 0 . 0 \ 

R == 0.0 0.0 0.0 (3.69) 

\ 0.0 0.0 0 . 0 / 

3.4.8 True of Date Ecliptic (TODEc) 

R/i and R/j for the TODEc system can be calculated 
using the following equations 

Rii = P r (i n )Rf (e)Ra (— A*F) (3.70) 

where t 0 is the epoch. Unlike the TOEEc sytem, for the 
TODEc system t 0 is a variable and usually comes from 
the epoch of the object whose state we wish to convert. 
P(to)is part of the FK5 reduction, algorithm and can be 
found in Vallado pgs. 214-219. e is given by Vallado, 1 
Eq. (3-52). A\F is given by Eq.(3-62) in Vallado. 1 For 
a more complete discussion, you can also refer to Seidel- 
mann, 2 pgs. 114 - 115. Because t 0 is not fixed for a 
MODEq system, the derivative of R/„; is non-zero. How- 
ever, we will assume its derivative is negligibly small so 
that 

/ 0.0 0.0 0 . 0 \ 

R := 0.0 0.0 0.0 (3.71) 

\0.0 0.0 0.0/ 

3.4.9 Celestial Body Fixed 

The body fixed coordinate system is referenced to the 
body equator and the prime meridian of the body. The 


body fixed system for Earth is found by using FK5 reduc- 
tion to the ITRF system as described by Vallado. The 
ITRF system is the earth fixed system. 

Vallado denotes the four rotation sequences required 
to transform from the ITRF to the FK5 system as [PM], 
the polar motion, [ST], the sidereal time, [NUT], the nu- 
tation, and [PR.ECj, the precession. GMAT calculates 
these four rotation matrices as described in Vallado. The 
rotation matrix from ITRF to FK5 can be written as fol- 
lows. 

R/i - P T N r ST T PM r (3.72) 

GMAT assumes that the the only intermediate rotation 
that has a significant time derivative is the sidereal time, 
[ST], So, we can write 

R/i = p t n t s't t pm t 

where ST is given by 

( ~u> e sin Oast a> e cos Oast 0.0 \ 

~(xj e cos Oast — w e sin Oast 0.0 J 

0.0 0.0 0 , 0 / 

and uj e is given by 

= 7.2921151467i)698e“ 5 ( 1 - L ° D - ) 

V 86400/ 

Note that the 2nd edition of Vallado 1 has inconsistencies 
in Eqs. (3.72) and (3.73) and they are discussed in the 
errata to the 2nd edition. We have modified equations 
Eqs. (3.72) and (3.73) according to the errata. 

For bodies other than the earth, the IAU gives the 
spin axis direction as a function of time with respect to 
the MJ2000Eq system and rotation of the prime merid- 
ian in the MJ2000Eq system. This data for all of the 
planets and many moons can be found in “Report of the 
IAU/IAG Working Group on Cartographic Coordinates 
and Rotational Elements of the Planets and Satellites: 
2000" Seidelmann 4 et.al. Figure 1 in this document ex- 
plains the three variables, a 0 , 0 o , and W , that are used 
to define the body spin axis and prime meridian location 
w/r/t J2000. The values of a Q , S 0 , and W for the nine 
planets and the Earth’s moon are found on pgs. 8 and 9. 

Using the notation found in the reference we can write 

Rh = Rs (90° + a 0 )Rf (90° - $ 0 )Rj (IV) (3.76) 


(3.73) 

(3.74) 

(3.75) 


For the derivative we assume that 


4 (R?(90° + a o )) = 
at 


and 


c it 


(R'f (90° - S 0 )) = 


/ 0.0 

0.0 

0.0\ 

0.0 

0.0 

0.0 

\0.0 

0.0 

0.0/ 

fo.o 

0.0 

0.0\ 

0.0 

0.0 

0.0 

lo.o 

0.0 

0.0/ 


(3.77) 


(3.78) 
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df 


(Rj(iV)) 


-tVsin(iV) 

-lVcos(!4'") 

0.0\ 

VVcos(lV’) 

— Vi'sin(H’) 

0.0 

0.0 

0.0 

0.0/ 


The remit is a rotation matrix, that is time invariant, for 
each celestial body. The individual matrices are shown 
below. 


(3.79) 

where W is the time derivative of W for the given body. 
Note, Seidelinann 4 does not provide the values for W so 
we include them in Table 3.1. In summary 


For the Sim 


R == R3 (90° + a 0 }Rf (90° 


3.4.10 Body Inertial 




(3.80) 


The Body Inertial axis system is an inertial system based 
on the equator of the celestial body chosen as the origin of 
the system. The origin of a Body Inertial system must 
be a celestial body, and cannot be a spacecraft, libration 
point etc. The axes are defined as follows (except for For Mercury 
Earth): 


Ru 

R 21 

«31 

Ru 

R-22 

R.% 2 

Ra 

R-23 

R 33 


• x-axis: Along the line formed by the intersection of 
the bodies equator and the x-y plane of the FK5 
system, at the J2000 epoch. 

• y- axis: Completes the right-handed set. 

• .z-axis: Along the bodies instantaneous spin axis 
direction at the J2000 epoch. 

For Earth, the Body Inertial axis system is the FK5 
system. For all other bodies, the Body Inertial axis sys- 
tem is based upon the bodies equator and spin axis at the 
J2000 epoch. So, Body Inertial is essentially a true-of- 
epoch system referenced to the chosen central body. The 
body orientation at the J2000 epoch is calculated from the 
IAU Data in Seidelinann 4 for the Sun, Mercury, Venus, 
Mars, Jupiter, Saturn, Uranus, Neptune, and Pluto. For 
the Moon, the orientation at the J2000 epoch comes from 
the DE405 files. Because the Body Inertial system is an 
inertial system, the derivative of the rotation matrix is 
always zero: 


For Venus 


Rfi = 



(3.81) 


R . (i 
R 21 
R 31 
Rl2 
R'22 
R 30 
RlS 

R-23 

R33 


R\i 

R-21 

R:n 

R\2 

R'22 

R32 

Rl3 

R23 

R33 


The rotation matrix, R M , is different for each celestial 
body. We begin by calculating the angles a and 6 used 
to define the bodies orientation with respect to the FK5 
system: 

a — °‘o{'lj20oo) (3.82) 

5 — 5 0 (Tj 2000) (3.83) 

Where T J200 0 - 2451544.9999999990686 TDB and the 
equations for a 0 and 5 0 are given by Seidelinann 4 and 
reproduced in Table 3.1. Finally, the rotation matrix is 
calculated using 


For the Earth 


0.9606338208497771 

0.2778176780544363 

0 

-0.2494239059807128 

0.8624542595398803 

0.4404093156676427 

0.1223534934723278 

-0.4230720836476433 

0.8977971010607901 


0.9815938660446788 

0.1909803187332696 

0 

-0.1677571842642237 

0.8622324234816705 

0.4779254910806334 

0.09127436261733374 

-0.4691287304711406 

0.8784003785150228 


0.9988399975085459 

0.04815245972043356 

0 

-0.04437694044018298 

0.9205233405740161 

0.3881590738545506 

0.01869081416890205 

-0.3877088083617988 

0.9215923900425705 


R J 3b ,i = R : ; (90° + a)Rj (90° - <5) 


(3.84) 


An - 
R 12 = 
RlS ~ 
R 2 I — 
R22 = 
R 23 — 

«31 = 

R 32 — 
R33 — 


(3.85) 


(3.86) 


(3.87) 


(3.88) 
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For the Moon 

R u = 0.998496505205088 
R 2 i = -0.0548154092680678 

#31 — 0 

R u = 0.0499357293985327 

R 22 - 0.909610125238044 

# 32 - 0.412451018902689 

Ru = -0.0226086714041825 

# 23 == ;: -0.411830900942613 

i ? 33 = 0.91097977859343 (3.89) 

For Mars 

Ru 
R21 
R31 
Ru 
R22 
R32 

#13 
R23 
#33 

For Jupiter 

R u = 0.9994209020316729 
#21 = -0.03402735050216735 

#31 - 0 

#i 2 = 0.0307100286015568 

R22 = 0.9019874904580493 

# 3 2 - 0.430668621100356 

#13 = -0.01465451212046692 

#23 = -0.4304192217768545 

# 33 - 0.9025101322420253 (3.91) 

For Saturn 

Ru - -0.6506284356468798 
R21 = 0.7593962330217962 
#31 — 0 

# 12 = -0.7545700815904118 

#22 = -0.6464935305479939 

#32 - 0.1125615694996715 

#13 = 0.08547883186107168 

#23 = 0.07323575787752883 

#33 = 0.9936447519469775 (3.92) 


For Uranus 

#11 = 0.9755767334083795 
#21 = -0.2196589111150187 

#31 — 0 

#12 = -0.05749969269512644 

#22 = -0.2553748540714755 

#32 - 0.9651308042166816 

#13 = -0.2119995815377986 

#23 = -0.9415591572895125 

#33 - -0.2617680858165513 (3.93) 

For Neptune 

#11 - 0.8717809455009272 
#21 = 0.4898958900230839 

#31 = 0 

# 12 = -0.3337976487639454 

#22 = 0.5940005535292547 

#32 = 0.7319443094160927 

#13 - 0.3585765089087282 

#23 = -0.6380951021167846 

#33 = 0.6813644604126335 (3.94) 

For Pluto 

#11 = 0.7311155947298647 
#21 - 0.6822536091093958 
#31 = 0 

#12 = -0.1077863335431382 

#22 = 0.1155058299435205 

#32 = 0.9874413955017208 

#13 = 0.6736854558650673 

#23 = -0.7219338031331282 

#33 = 0.1579857286263988 (3.95) 

3.4.11 Object Referenced 

An object referenced system is a coordinate system whose 
axes are defined by the motion of one object with respect 
to another object. GMAT allows the user to define many 
different types of Object Referenced system. In Fig. 3.6 
we see a diagram that defines the directions a user can 
choose from in creating an Object Referenced coordinate 
system. There are six directions. One is the relative po- 
sition, denoted here by r, of the secondary object with 
respect to the primary object, expressed in an inertial 
frame. The second is the the relative velocity, denoted 
here by v, of the secondary object with respect to the 
primary, expressed in an inertial frame. The third direc- 
tion is the vector normal to the direction of motion which 


= 0.6732521982472343 
= 0.7394129276360177 
= 0 

= -0.5896387605430038 
= 0.5368794307891334 
= 0.6033958972853944 
= 0.4461587269353554 
= -0.4062376142607542 
= 0.7974417791532832 (3.90) 
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Figure 3.6: Diagram of an Object Referenced Coordinate System 


is denoted by n and is calculated from n = r x v. The 
remaining three directions are the negative of of the first 
three. 

In GMAT, a user can use the directions described 
above to define an Object Referenced coordinate system. 
In doing so, the user can choose two of the available di- 
rections. and define which of the three axes, x, y, and z, 
they desire the direction to be aligned with. Given this 
information, GMAT automatically constructs an orthog- 
onal coordinate system. For example, if user chooses the 
x-axis to be in the direction of r and the 2 -axis to be in 
the direction of n, the GMAT completes the right-handed 
set by setting the y-axis in the direction of n x r. Obvi- 
ously there are some choices that not. yield an orthogonal 
system, or that yield a left handed system. GMAT does 
not allow the user to select these pairs of axes and throws 


and 


xxytxxy 


(3.99) 


If the user defines the y-axis and 2 -axis, then GMAT de- 
termines the z axis using 

x -yxz (3.100) 

and 

x = yxz + yxz (3.101) 

And finally, if the user defines the x-axis and 2 -axis then 
GMAT determines the z axis using 

y = z x x (3.102) 

and 

y — zxx + zxx (3.1 03) 


an error message. 

In general, given the unit vectors that define the axes 
system of the Object Referenced system, but expressed 
in the inertial frame, GMAT uses the following equations 
to determine and R/, . 

R(i = [ x y z ] (3.96) 

R/,: = [ * 9 ft ] (3.97) 

Recall that the user chooses two axes to an Object 
Referenced system among the following choices: r , v, ft, 
— r, — v, and — n. In general, one of the axes chosen by 
the user must be either h, or — n. 

If the user defines the a:- ax is and y-axis then GMAT 
determines the z axis using 

z — x. x y (3.98) 


Depending on the users choice of axes for an Object 
Referenced coordinate system, GMAT will need to calcu- 
late r, v, h, r, v, and n. These are given by: 


r = 

r r 

IMI ~ r 

(3.104) 

v — 

V V 

IMI ~ v 

(3.105) 

n = 

r x v 

(3.106) 

! r x v |j 

Jf v 

r = - 
r 

- - (r ■ v) 
r 

(3.107) 

• a 
v — — • 

V 

V . 
tv • ai 

V 

(3.108) 

r x a 

n , _ 

f. r x a • n) 

71 

(3.109) 

n 
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3.4.12 Geocentric Solar Ecliptic (GSE) 

The Geocentric Solar Ecliptic system is a time varying 
axis system often used to describe and analyze the Earth’s 
magnetic field.. The coordinate system is defined such 
that 

Tsun ( 3 . 110 ) 


X : 


where r sura is the vector from the Earth to the Sun in the 
MJ20C)0Eq axis system. The 2 -axis is defined to be the 
ecliptic pole. To ensure we have an orthogonal system, 
we calculate z using 


x v*. 


P'sitn ^ 


(3.111) 


(3.119) 

(j> d = 78.6° N (3.120) 

The dipole vector in the Earth Fixed system is simply 

( cos 4>d cos (— A d )\ 
cos 4>d sin (—Ad) (3.121) 

sin 4> d j 

If Rj F is the rotation matrix from the Earth Fixed 
frame to MJ2000Eq at the current epoch, then we can 
write the vector that describes the dipole direction in the 
inertial frame as 


Finally, the y-axis completes the right-handed set 

y — z x x (3.112) 

We can construct the rotation matrix that goes from the 
GSE axis system to the MJ2000Eq axis system as 


{r rf }, = R if { r d } r 
Then, the y-axis is defined as 

- __ {r d }j x x 


{ r d}i x x|| 


R.;,: = 


(3.113) 


the z-axis is defined as 


We also need to compute the derivative of the rotation 
matrix. We start by computing 


dx __ v sur 
di r su „ 


(3-114) 


and 


z — x. x y 


R /» = [ x y z ] 


13.122) 


(3.123) 


(3.124) 


(3.125) 


where v stm . is the velocity of the Sim with respect to the 
Earth in the MJ2000Eq system. We can approximate the 
derivative of the z axis using 


To calculate the derivative of the rotation matrix, we 
know that 


dz 

dt 


0 


R a 


dy „ dx 
-- = z x — 
dt dt 

die dy dz 

dt dt dt 


(3.115) 

(3.116) 

(3.117) 


dx 


dt r. 


f 




Let’s define 


and 


y = (Riv{rd}pj x x 

V= II (^/f{ r <*}*=') x xj 


(3.126) 

(3.127) 

(3.128) 


then 


3.4.13 Geocentric Solar Magnetic (GSM) 


(3.118) 


dy 

dt 


= y= (Hi pil'd} f J X x + (Rjf{ r d } F ) x (3.129) 
Now we can write 


Let’s define the spherical coordinates of the Earth’s 
dipole in the Earth fixed frame to be A d and <p d - The 
location of the dipole actually changes with time. Also, 
the dipole does not actually pass through the center of 
the Earth. However, GMAT currently assumes that the 
dipole direction is constant, and passes directly through 
the center of the Earth. If this approximation is not suffi- 
cient for future studies, the model will have to be updated. 


Finally, 


and 


dy i y „ /„ y 

Vi =y= — y y • - 

dt y \ y 


z=xxy+xxy 


= 


dx dy dz 

dt dt dt 


(3.130) 


(3.131) 


(3.132) 



3.5. APPENDIX 1: DERIVATIVES OF OBJECTREFERENCED UNIT VECTORS 


29 


3.5 Appendix 1: Derivatives of Ob- 
jectReferenced Unit Vectors 

The derivations of the above quantities are shown below. 
We start by deriving two derivatives with respect to n, 
where n is given by: 

n = r x v (3.133) 

We need to determine two derivatives of n. First 

(3.134) 


• dv d , , a v /v • a\ 

v “S-S (vl 'T = ) 

which can be rewritten as 


a a v 
v — — — — (v ■ a) 


Finally, 


x « / _ i \ r x a rt , . 

11 = a7 ( nn ) = W^-^( rxa ' n ) 


dn 


d , .dr dv 

¥ = S (rxv > = S xv+r>< W 


dt 


n n- 




x r x a 

11 / ~ N 

n = 

— (r x a • n) 

n 

n 


so we know that 


dn 

dt 


— r x a 


The next useful derivative is 

dn d|| nil d , T . . 1/2 n T dn 

— = — — — — — inn) — 

dt dt dt' n dt 

So we can write 


dn n . . 

— = - • (r x a) 
dt n 


The following two derivatives are also useful 
diirll 


dr 

dt 


■■ ■■ = —(rT r \ 1/2 _ 111 

dt dt ' 1 r 


dv d 
dt c 
so we can write 


dr v • r 
dt r 

d 
dt 


- -ixCUvU 2 - 


v • a 


dv v ■ a 
dt v 


(3.135) 

(3.136) 

(3.137) 

(3.138) 

(3.139) 

(3.140) 

(3.141) 


(3.147) 


(3.148) 


(3.149) 

(3.150) 


v 



r 


„ r x v 

n = 

IT x v|| 

The time derivatives are derived as follows. 



which can be rewritten as 


dr 

dt 


v r ., 
(r • v) 

r r 


(3.142) 

(3.143) 

(3.144) 

(3.145) 


(3.146) 
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Tabl e 3.1: Recommended Values for Pole and Prime Meridian Locations of the Sim and Plan ets 4 

Name Values 

Sun a Q — 286.13° (deg) 

S 0 = 63.87° (deg) 

W ------ 84.10° -1- 14.1844000°flf (deg) 

W = 14.1844000° (deg/s) 


Mercury a a — 281.01 — 0.0337' 

S a — 61.45 — 0.0057’ 

W = 329.548 + 6.1385025d 
W = 6.1385025 


Venus a 0 = 272.76 

S a = 67.16 

W = 160.20 - 1.481 3688(7 
W = -1.4813688 

Earth a 0 — 0.00 — 0.6417' 

5 0 = 90.00 - 0.5577’ 

W = 190.147+360.9856235 d 
W = 360.9856235 

Earth Data is included for completeness only, GMAT uses FK5 reduction 
for the Earth 

Mars a 0 = 317.68143 - 0.10617’ 

<5„ = 52.88650 - 0.06097’ 

IV : 176.630 + 350. 89198226c/ 

W = 350.89198226 

Jupiter a„ ----- 268.05 — 0.0097' 

S 0 — 64.49 + 0.0037’ 

VV = 284.95 + 870.5366420d 
W ==: 870.5366420 

Saturn a 0 = 40.589 — 0.0367’ 

S 0 = 83.537- 0.0047' 

W = 38.90 + 810.7939024 d 
W = 810.7939024 

Uranus a Q — 257.311 
6 0 = -15.175 

W = 203.81 - 501.1600928d 
W m -501.1600928 

Neptune a n = 299.36 + 0.70 sin TV 
S 0 ------ 43.46 — 0.51 cos TV 

W = 253.18 + 536.3128492(7 - 0.48 sin N 
W = 536.3128492 - 0.48TV cos N 
TV = 357.85 + 52.3167’ 

TV = 6.0551 x 10 -4 (deg/day) 

Pluto a c , — 313.02 
5 0 = 9.09 

W — 236.77 — 56.3623195(7 
W = -56.3623195 
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Table 3.2: Recommended Values for Pole and Prime Meridian Locations of Luna 4 

Name Values 

Luna 


- 269.9949 +0.00317' -3.8787 sinAT -0.1 204 sin £2 

+0.0700 sin A3 -0.0172 sin E4 +0.0072 sin E6 

—0.0052 sin A10 —0.0043 sinAT 3 


S 0 = 66.5392 


+0.01307’ 
-0.0278 cos A3 
+0.0009 cos A'7 


+1.5419 cos El 
+0.0068 cos A'4 
+.0008 cos AT 0 


+0.0239 cos E2 
-0.00292 cos.A'6 
-0.0009 cos AT 3 


where 


W = 38.3213 


+13.17635815(7 
+0.1208 sin A2 
+0. 0252 sin Eh 
-0.0046 sin A8 
+0.0040 sin Al l 


-1.4 x 10 i2 e/ 2 
-0.0642 sin A3 
—0.0066 sin A6 
+0.0028 sin A9 
+0.0019 sin AT 2 


+3.5610 sin AT 
+0.0158 sin A'4 
—0.0047 sin A7 
+0.0052 sin AT 0 
-0.0044 sin AT 3 


W = 


+13.17635815 
-.01280 cos A2 
+.0248 cos A'5 
-.0015 cos A8 
+.00001 cos All 


-2.8 x 10 _12 d 
—.835 cos A3 
—.17 cosA6 
+.0049 cos A9 
+.00031 cos AT2 


-.18870 cos AT 
+.211 cos A4 
—.061 cos A7 
— .00083 cos A10 
— .057 cos AT 3 


A1 = 125.045 - 0.0529921c/ 
A3 = 260.008+ 13.0120009(7 
A5 - 357.529 + 0.9856003a 
El - 134.963 + 13.0649930(7 
A9 = 34.226+ 1.7484877(7 
All - 119.743 + 0.0036096(7 
AT 3 - 25.053 + 12.9590088(7 


A2 = 250.089 - 0.1059842(7 
A'4 = 176.625 + 13.3407154(7 
A'6 - 311.589 + 26.4057084(7 
A8 - 276.617 + 0.3287146(7 
AlO = 15.134 - 0.1589763(7 
A12 — 239.961 + 0.1643573(7 
















Chapter 4 

Calculation Objects 


GMAT has the ability to calculate numerous quanti- 
ties that are dependent upon the states of objects, co- 
ordinate systems, and the mission sequence. These cal- 
culation objects can range from the spacecraft state, to 
the local atmospheric density, to the positions of celestial 
bodies with respect to spacecraft, or other celestial bod- 
ies. In chapter, we present how GMAT performs these 
calculations by showing the mathematical algorithms. 

The chapter begins by discussing different orbit state 
representations. Each of the orbit state representations 
available in GMAT are defined. Next we present the algo- 
rithms used to convert between different state represen- 
tations. These include the Keplerian elements, modified 
Keplerian elements, Cartesian state, spherical state, and 
the equinoctial elements. In the second section we present 
how GMAT calculates all calculation parameters. Exam- 
ples include the orbit period, percent shadow, and energy. 
The algorithms to calculate all parameters are included 
and described in detail. We conclude this chapter with a 
presentation of the algorithms used to calculate libration 
point and baxycenter locations. 

4.1 Spacecraft State Representa- 
tions 

There are several state representations that can be used 
in GMAT to define the state of a spacecraft object. These 
include the Keplerian elements, Cartesian state, equinoc- 
tial elements, spherical elements, and the modified Kep- 
lerian elements. In the following subsections, we discuss 
the definitions of these states types, and show how GMAT 
converts between the different state representations. 


4.1.1 Definitions 

The Keplerian elements are one of the most commonly 
used state representations. They provide a way to de- 
fine the spacecraft state in way that provides an intuitive 
understanding of the motion of spacecraft in orbit. The 


Keplerian elements are denoted a, e, i, w, fi, and v. They 
are defined in detail in Table 4.2 and illustrated in Fig. 
4.1. Sections 4.1.2 and 4.1.3 show the algorithms that 
GMAT uses to convert between the Keplerian elements 
and the cartesian state. 

The cartesian state is another common state repre- 
sentation and is often used in the numerical integration 
of the equations of motion. The cartesian state with re- 
spect. to a given coordinate system is described in detail 
in Table 4.1. 

The equinoctial elements axe a set of non-singular el- 
ements that can be used to describe the state of a space- 
craft. Because they are nonsingular, they are useful for 
expressing the equations of motion in Variation of Pa- 
rameters (VOP) form. The elements can be unintuitive 
to use however. The equinoctial elements are described 
in detail in Table 4.4. 

The modified Keplerian elements are similar to the 
Keplerian elements except a and e are replaced with the 
radius of periapsis r p , and the radius of apoapsis r a , r p 
and r a are often more convenient and intuitive for describ- 
ing the dimensions of a Keplerian orbit than a and e. The 
modified Keplerian elements are defined in detail in Table 
4.6. Note that both the Keplerian and modified Keple- 
rian elements are undefined for parabolic orbits because 
the semimajor axis is infinite. Currently, GMAT does not 
support parabolic orbits for this reason. Let’s begin by- 
looking at how GMAT converts from the Cartesian state 
to the Keplerian elements. 


4.1.2 Cartesian State to Keplerian Ele- 
ments 

The conversion from the Cartesian state to the Keplerian 
elements has four special cases: elliptic inclined, circular 
inclined, elliptic equatorial, and circular equatorial. Cer- 
tain orbital elements are undefined for some of the cases. 
For example, the right ascension of the ascending node, 
fi, is undefined for equatorial orbits. However, computer 
systems don’t handle undefined values gracefully. In this 
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Figure 4.1: The Keplerian Elements 

section, we’ll see how the orbital elements are defined for Recall that for parabolic orbits, the semimajor axis is 
each of the special cases, and how GMAT calculates the infinite and the energy is zero. GMAT checks to see if 
orbital elements for each case. the orbit is near parabolic and returns an error message 

in this case. The following error check avoids the possi- 
Given. r, v , and p bility of a divide by zero: 


Find: a, e, i, w, ft, and v 

We begin by calculating the specific angular momen- 
tum and its magnitude. 


if |1 -e\> 1CT 30 , then 


« = 

2 € 


h = r x v 


h = El h 


(4.1) otherwise, report error and return. The error reported is: 
“Warning: GMAT does not support parabolic orbits in 

(4.2) conversion from cartesian to Keplerian state” . 


A vector pointing' in the direction of the line of nodes is If the orbit is not parabolic according to the above 

defini tion, then we continue and calculate the inclination. 

n = [ 0 0 1 j T x h (4.3) 


The orbit eccentricity and energy are calculated using 

There are four special cases for ft, u>, and v and each 
r _ |j r || (4 5} case is treated differently. 


(v 2 — --)r — (r ■ v)v 


V" fj. 


. , . Special Case 1: Non-circular, Inclined Orbit 
(4.6) 

if (e > 10 _u ) and ( i > 10 -11 ), then 

( 4 -7) n = co8- 1 ^) 

\ V< / 

(^•8) Fix quadrant for ft: if n y < 0, then ft — 2 tt -- ft 
(4.9) u> == cos" 1 

V ne / 
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Table 4.1: The Cartesian State 


Symbol 

Description 

X 

x-component of position 

y 

y-eomponent of position 

Z 

2 -component of position 

X 

x-component of velocity 

y 

y-component of velocity 

Z 

z-component of velocity 


Table 4.2: The Keplerian Elements (also see Fig. 4.1) 


Symbol 

Name 

Description 

a 

semimajor axis 

The semimajor contains information on the type and size 
of an orbit. If a > 0 the orbit is elliptic. If a < 0 the orbit 
is hyperbolic, a — co for parabolic orbits. 

e 

eccentricity 

The eccentricity contains information on the shape of an 
orbit. If e = 0, then the orbit is circular. If 0 < e < 1 the 
orbit- is elliptical. If e — 1 the orbit is parabolic. If e > 1 
then the orbit is hyperbolic. 

i 

inclination 

The inclination is the angle between the 2/ axis and the 
orbit normal direction h. If i < 90° then the orbit is 
prograde. If i > 90° then the orbit is retrograde. 

u? 

argument of periapsis 

The argument of periapsis is the angle between a vector 
pointing at periapsis and a vector pointing in the direction 
of the line of nodes. The argument of periapsis is undefined 
for circular orbits. 

n 

right ascension of the as- ft is defined as the angle between xj and N measured coun- 
cending node terclockwise. N is defined as the vector pointing from 

the center of the central body to the spacecraft, when 
the spacecraft crosses the bodies equatorial plane from the 
southern to the northern hemisphere. Q is undefined for 
equatorial orbits. 

V 

true anomaly 

The true anomaly is defined as the angle between a vector 
pointing at periapsis and a vector pointing at the space- 
craft. The true anomaly is undefined for circular orbits. 


Table 4.3: 

Keplerian Elements for Special Cases 


Orbit Type 

Numerical Threshold 

Description 

Elliptic Inclined 

e > IQ" 11 , i > 10” u 

Q is the angle between the x-axis and the fine of 
nodes, w is the angle between the line of nodes and 
the eccentricity vector, u is the angle between the ec- 
centricity vector and the spacecraft position vector. 

Elliptic Equatorial 

e > 1(T S! , i < 10 -11 

ft = 0, cn is the angle between the x-axis and the 
eccentricity vector, i> is the angles between the eccen- 
tricity vector and the spacecraft position vector. 

Circular Inclined 

e < 10~ n . i > 10 _u 

ft is the angle between the x-axis and the fine of 
nodes, w — 0, v is the; angle between the line of nodes 
and the spacecraft position vector. 

Circular Equatorial 

e < IQ -11 , i < IQ- 11 

ft = 0, w = 0, v is the angle between the x-axis and 
the spacecraft position vector. 
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Fix quadrant for w: if e z < 0, then u> — 2?r — ui 



Fix quadrant for v: if r • v < 0, then v — ‘hx — v 
Special Case 2: Non-circular, Equatorial Orbit 
if (e > 10” n ) and ( i < 10" 11 ), then 

n = o 


(4.14) 


if (e < 10 n ) and (i < 10 11 ), then 


Q = 0 

LJ = 0 



Fix quadrant for u: if r y < 0, then v — 2r; — v 


(4.21) 

(4.22) 

(4.23) 


In the next section, we look at how to perform the 
inverse transformation and convert from Keplerian ele- 
^ ’ J merits to the Cartesian state vector. 


i0 = cos' 1 — 14.16) 

e 

Fix quadrant for w: if e y < 0, then w = 2vr — w 

u cos -1 (4.17) 

V er J 

Fix quadrant for v: if r ■ v < 0, then i/ = — u 

Special Case 3: Circular , Inclined, Orbit 
if (e < 10 " 11 ) and ( i > 10"' 11 ), then 

n = cos _1 (~') (4.18) 

\ n / 

Fix quadrant, for fl: if n y < 0, then fi = 2rr — ft 


4.1.3 Keplerian Elements to Cartesian 
State 

The transformation from the Keplerian elements to the 
cartesian state is one of the most common state transfor- 
mations in astrodynamics. We previously defined both 
state; types and refer you to Tables 4.1 and 4.2 for their 
definitions. Below we show the algorithm that GMAT 
uses to convert from the Keplerian elements to the Carte- 
sian state. 

Give: a, e, i, Q, cu, //, and p 
Find: r and v 


uj = 0 



Fix quadrant for v: if r z < 0, then v = 2tt — v 
Special Case Circular, Equatorial Orbit 


(4.19) We begin by checking that the magnitude of the po- 
sition vector is not infinite. This avoids the possibility of 

(4.20) a division by zero. 

if (1 + ecos v < le-30), then display error and return: 
“Warning: Radius is near infinite in keplerian to cartesian 
conversion” 
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If their is not a divide by zero issue we continue by cal- 
culating the semilatus rectum, and the radius. 

p = a(l - e 2 ); (4.24) 


We begin by using the mean longitude, A, to find the 
true longitude F. The equation relating the two is tran- 
scendental: 


A = F + h cos F — k sin F (4. 32) 


r - ~ (4.25) 

1 + e cos v 

The position components of the cartesian state vector are 
calculated using the following three equations. 

x — r (cos (u! + v ) cos ft — cos i sin (u; + is) sin ft) (4.26) 

y — r (cos (u + is) sin ft 4- cos i sin (w + is) cos ft) (4.27) 
z = r (sin (a) + is) sini) (4.28) 


We use the Newton-Raphson method to solve for F, using 
A as the initial guess. We iterate on the following equation 
until \F(i + 1) - F{i ) j < IQ- 10 . 

F(t + 1) = F(i) - (4.33) 

where 

f(F) = F + hcos(F) — ksin(F) — A (4.34) 

f'(F) == 1 — h. sin(i'’) — k cos{F) (4.35) 


Before calculating the velocity components we check 
to ensure the orbit is not parabolic. This avoids another 
possible division by zero. 

if (INI < le — 30), then error and return: “Warning: 
GMAT does not support parabolic orbits in conversion 
from keplerian to cartesian elements" . 

If the orbit is not parabolic, we continue and calculate 
the velocity components using 


I n 

x — if — (cos is + e) (— sin w cos ft — cos i sin. ft cos u>) - 

V V ' 


,/^"sin n(cosu,'Cosft — cos* sin ft since?) 

V p 


(4.29) 


/ .. 

V — \l (cos v + e ) (~ sin lj sin ft + cos i cos ft cos to) ■ 


V P 


I u 

i I --- sin v (cos u sin ft -f cos i cos ft sin u ) 

V P 


(4.30) 


: = ^j— |(cos;z + e) sin i cos w — sin v sini sin w] (4.31 


Now let’s look at how to calculate the cartesian state 
given the equinoctial elements. 


4.1.4 Equinoctial Elements to Cartesian 
State 

The equinoctial elements used in GMAT are defined in 
Table 4.4. The algorithm to convert from equinoctial el- 
ements to the cartesian state was taken from the GTDS 
Mathematical Theory . 5 

Given: a, h , k. p, q, A, and u 

Find: r and v 


Once the true longitude is calculated, we continue 
with 


»- —L 

1 + -s/1 - k 2 - k 2 
mu 

n = 3 " 

a 6 

(4.36) 

(4.37) 

r = a(l — k cos F — h sin /) 

(4.38) 

The cartesian components expressed in the 
coordinate system can be calculated using. 

equinoctial 

X\ --a [(1 — h 2 0) cos F + hk(i sin F - k ] 

(4.39) 

Y\ — a [(1 — k 2 0) sin F + hk3 cos F — h\ 

(4.40) 

X, - n °T \hki3 cos F (1 H 2 j3) sin F) 

r ' - 

(4.41) 

Vi - ^ [(1 - k 2 0) cos F - hkp sin F] 

(4.42) 

The transformation from the equinoctial system to the 
inertial Cartesian system is given by 

r - Xj + Fig 

(4.43) 

v = ATf + Ng 

(4.44) 

rvhere 

[ f * W ]-l + p?. + ,2 Q 

(4.45) 


and 


/ 1 -- p 2 -f q 2 2pqj 2 p 

Q = 2 pq (1 +p 2 -q 2 )j -2 q 

\ - 2 pi 2 q (l -p 2 ~q 2 )j 

(4.46) 

and j = 1 for direct orbits ( 0 < i < 90° ) 
j = -1 for retrograde orbits ( 90 < i < 180° ) 


Now let’s look at how to calculate the cartesian state 
given the equinoctial elements. 
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Table 4.4: The Equinoctial Elements 


Symbol 

Description 


a 

The semimajor contains information on the type and 
orbit, is elliptic. If a < 0 the orbit is hyperbolic. 

size of an orbit. If a > 0 the 

h 

The projection of the eccentricity vector onto the y ep 

axis. 

k 

The projection of the eccentricity vector onto the it ep 

axis. 

V 

The projection of N onto the y ep axis. 


q 

The projection of N onto the Xq> axis. 


A 

The mean longitude. 



4.1.5 Cartesian State to Equinoctial 
Elements 

The equinoctial elements used in GMAT are defined in 
Table 4.4. The algorithm to convert from the carte- 
sian state to the equinoctial elements was taken from the 
GTDS Mathematical Theory. 5 

Given: r, v, and ft 


k = e • f 

(4.57) 

kx 

P ------ , A 

1 + hi 

(4.58) 

Hqj 

q= -l + i4 

(4.59) 


The final element to calculate is the mean longitude, 
A. We begin by computing the true longitude, F, using 


Fine: a, h, k, p , q, A, and u 

We begin by calculating the semimajor axis and the 
eceentriei ty vector. 

r = |jr[| (4.47) 


v = ll v ll 

(4.48) 

1 

Cl •■■■ n 

2 v 2 

(4.49) 

r ji 


li 

! 

-1 \ *-! 

\ 

'►T 

X 

X 

< 

(4.50) 

The angular momentum vector is 


- r x v 

|ir x v|| 

(4.51) 

The unit vectors that define the equinoctial coordinate 
system can be calculated using 

Ir* 

ii 

i— 1 

i 

* 

N '-J. 

(4.52) 

i 

ii 

(4.53) 

!, = -hi 

(4.54) 

where j is described in Section 4.1.4. 


g = h x f 

(4.55) 


We now have the necessary information to calculate 
the elements h, k, p, and q using the following relation- 
ships. 

(4.56) 


and 


cos F — k + 
sin F - h + 


Xl T • f 

(4.60) 

*i = r-g 

(4.61) 

(1 - k 2 [3) Xl - hkpYi 

ay / 1 — h 2 — fc 2 

(4.62) 

(1 - F 2 0) Y X - KkpXy 

(4.63) 

ay 1 — h? — k 2 

. /sin F\ 

-*"M CCS r) 

(4.64) 


where 0 is given by Eq. 4.36. The mean longitude is 
computed using the generalized Kepler equation 


A — F + h cos f — k sin F 


(4.65) 


Nov/ let’s look at transformations involving the spher- 
ical elements. 


4.1.6 Cartesian State to Spherical AZFPA 
State 

The spherical state, with azimuth, a/, and flight path 
angle, ?/>, is described in Table 4.5 and Fig. 4.2. The al- 
gorithm below shows how GMAT converts from the carte- 
sian state to the spherical state with azimuth and flight 
path angle. 

Given: r and v 

Find: r, A, <5, v, ijj, and aj 


h = e- g 
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Table 4.5: The Spherical Elements 


Symbol 

Name 

Description 

r 

V 

Magnitude of the position vector, ||rij 

A 

Right Ascension 

The between the projection of r into the x — y and the 
x-axis. Measured counterclockwise. 

5 

Declination 

The angle between r and the x — y plane measured in the 
plane formed by r and z . 

V 

V 

Magnitude of the velocity vector, ||v|| 

y; 

Vertical flight path angle 

The angle measured from a plane normal to r to the; ve- 
locity vector v, measured in the plane formed by r and 

ay 

Flight path azimuth 

The angle measured from vector perpendicular r and 
pointing north, to the projection of v into a plane nor- 
mal to r. 


We begin by calculating the right ascension A, and the 
declination 6. 

r = |jr|| (4.66) 

A = tan 2 : (y, x ) (4.67) 

d = sin“ 1 (-) (4.68) 

r 

The magnitude of the velocity vector is simply 

v = 1 1 v 1 1 (4.69) 

We calculate the vertical flight path angle, psi, using 

ip = cos" 1 f— ) 14.70) 

\ TV / 

To calculate the azimuth angle, a,, we begin by calcu- 
lating the rotation matrix from the frame in which the 
cartesian state is expressed in, T, to a local frame, Ti, 
where z is a unit vector that points north. The basis 
vectors of Tt expressed in Ti can be calculated using 


Now that we have looked at how to convert from the 
Cartesian state to the spherical state, let’s look at the 
inverse transformation that converts from the spherical 
state (with ip and ay) to the cartesian state. 


4.1.7 SphericalAZFPA State to Cartesian 
State 

In this section we present the algorithm used to convert 
from the spherical state (with ip and ay) to the cartesian 
state. 

Given: r, A, 5, v, ip, and ay 
Find: r and v 

The components of the position vector are calculated 
using 


/cos(d) cos(A)\ 



x = r cos 5 cos A 

(4.77) 

[ cos(d) sin(A) j 

(4.71) 


y = r cos d sin A 

(4.78) 

V sin(d) j 



z — r sin d 

(4.79) 

/ cos(A 7r/2)\ 





[ sin(A + 7 t/2) 

(4.72) 

We 

can write the velocity vector in terms 

of v, ’ip, and 

V o J 


ctf as. 



sin(d) cos(A)' 1 

\ 

v — 

v fcos(y;)x + sin(y6) sin(oy)y + sin (ip) t 

:os(ay)zj 

— sin(d) sin.(A) 

(4.73) 



(4.80) 

„ cos(d>) , 

/ 

where, 

x, y, and x are given in Eqs. (4.71), 

(4.72), and 


We can write the tranformation matrix that goes from Ti 
to Tt, Rri, as 

= [ x y z ] T (4-74) 

The velocity in the local frame, v 1 , can be written as 

v' = R«v (4.75) 

Finally, we calculate the azimuth angle using 

ay — tan7 x (v' y ,v' z ) (4-76) 


(4.73) respectively. Breaking down Eq. (4.80) into com- 
ponents gives us 


v x — v [cos ii> cos (I cos A— sin ip (sin a j sin A 
+ cos ay sin d cos A)] 

v y = vicos ip cos 5 sin A+ sin </>(sin a j cos A 


(4.81) 

• x ■ w ( 4 - 82 ) 
— cos ay sm o sin A J : 

v z — t>[cos ip sin d + sin ip cos arcosS ] (4.83) 
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4.1.8 Cartesian State to SphericalRADEC 
State 

The conversion form the Cartesian state to the spherical 
state with right ascension of velocity, A„ , and declination 
of velocity, S v , is very similar to the transformation shown 
in Sec. 4.1.6. The algorithm to calculate A.„. and 5 V is 
shown below. 

Given: r and v 

Find: r, A, 5, v, A„, and 6 V 

To calculate r, A, 5 , and v we use Eqs. (4.66), (4.67), 
(4.68), and (4.69) respectively. The right ascension of 
velocity, A„, and declination of veloci ty, S v . are calculated 
using 

At, = taxi 2 1 (v y , v x ) (4.84) 

5 V = sin -1 !"--) (4.85) 

v 

In the next section, we show the the transformation 
from the spherical state with right ascension of velocity, 
A„, and declination of velocity, S v , to the cartesian state. 

4.1.9 SphericalRADEC State to Cartesian 
State 

This transformation is similar to the conversion presented 
in Sec 4.1.7. The primary difference is how the velocity 
is represented. 

Given: r, A, S, v, A„, and $ v 

Find: r and v 

The position components are calculated using Eqs. 


(4.77), (4.78), 

and (4.79). 

The velocity 

components are 

calculated using 




v x ------ 

v cos A cos 5 

(4.86} 


Vy - 

v x tan A 

(4.87) 


v z = 

v sin 5 

(4.88) 


the radius of periapsis, r p , to describe the size and shape 
of an orbit. The remaining elements, i, fi, w, and v, are 
the same for both the Keplerian and modified Keplerian 
elements. The modified Keplerian elements, like the Ke- 
plerian elements, are undefined for parabolic orbits. Let’s 
look at how GMAT calculates the modified 

Given: a. e, i, u>, fi, and or r, v, and /( 

Find: r p and r a 

If we are given the Cartesian state, we first calculate 
the orbital elements using the algorithm in Sec. 4.1.3. 
Knowing the Keplerian elements, we calculate r 0 and r P 
using 

r a = a(l + e) (4.89) 

r. p = a(l — e) (4.90) 

Now let’s look at the inverse transformation. 


4.1.11 Modified Keplerian Elements to Ke- 
plerian Elements 

The conversion from modified Keplerian elements to the 
Keplerian elements is discussed below. To perform the 
conversion, we use relationships that allow us to right the 
semimajor axis, a, and the eccentricity, e, in terms of r a 
and r p . 

Given: r v , r a , i, u>, fi, and v 
Find: a and e 

We begin by calculating the eccentricity using 


The semimajor axis is calculated using 


In the last few subsections, we have looked at trans- 
formations involving the spherical elements. Now let’s 
look at transformations involving the modified Keplerian 
elements. 

4.1.10 Keplerian or Cartesian, to Modi- 
fied Keplerian Elements 

The modified Keplerian elements, described in Table 4.6, 
are similar to the classical Keplerian elements. The modi- 
fied Keplerian elements use the radius of apoapsis, r a , and 



This concludes our discussion of state transformations. 
In the last few subsections we presented the algorithms 
used to convert between different orbit state representa- 
tions used in GMAT. These include the Cartesian state, 
the Keplerian elements, the modified Keplerian elements, 
and two spherical state parameterizations. In the next 
section, we present the algorithms used to calculate prop- 
erties such as orbit period, beta angle, and mean motion 
to name a few. 
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Table 4.6: The Modified Keplerian Elements 


Symbol 

Name 

Description 

r P 

radius of periapsis 

The radius of periapsis is the radius at the spacecrafts clos- 
est approach to the central body. The radius of periapsis 
must be greater than zero, parabolic orbits are not cur- 
rently supported. 

t'a 

radius of apo apsis 

For an elliptic orbit r a is the radius at the spacecrafts far- 
thest distance from the central body and r Q > r v . For 
hyperbolic orbits, r a < r p and r a < 0 

i 

inclination 

The inclination is the angle between the : z; axis and the 
orbit normal direction h. If i < 90° then the orbit is 
prograde. If i > 90° then the orbit is retrograde. 

UJ 

argument of periapsis 

The argument of periapsis is the angle between a vector 
pointing at periapsis, Xj,, and a vector pointing at the 
spacecraft. The argument of periapsis is undefined for cir- 
cular orbits. 

n 

right ascension of the as- 
cending node 

Q is defined as the angle between x/ and N measured coun- 
terclockwise. N is defined as the vector pointing from 
the center of the central body to the spacecraft, when 
the spacecraft crosses the bodies equatorial plane from the 
southern to the northern hemisphere. Q is undefined for 
equatorial orbits. 

V 

true anomaly 

The true anomaly is defined as the angle between a vec- 
tor pointing at periapsis, x p , and a vector pointing at the 
spacecraft . The true anomaly is undefined for circular or- 
bits. 


4.2 Simple Parameters 


Simple parameters, which we will abbreviate as simply 
“parameters” , are properties of spacecraft or other ob- 
jects that are only dependent upon one of the following: 
CoordinateSystem, CentralBody, or None. An example 
of a simple parameter is the magnitude of a spacecrafts 
velocity vector. The spacecrafts velocity vector is depen- 
dent upon the coordinate system in which it is expressed. 
Once we have specified a coordinate system, it is trivial to 
calculate the velocity vector, and therefore its magnitude, 
in that coordinate system. 

In GMAT, the syntax to specify a simple parameter 
is 

Obj ectName . Dependency . ParameterName 

So, to calculate the magnitude of the velocity, of a space- 
craft named Sat, in the Earth Fixed frame, we would use 

Sat . EarthF ixed . VMAG 

GMAT has the ability to calculate many parameters in 
addition to VMAG. In the following subsections, we present 
the algorithms used to calculate all parameters in GMAT. 
We begin each subsection with a description of the param- 
eter, and then give the type of dependency. 


4.2.1 AlGregorian 

Description: AlGregorian is the epoch of an object, in 
the Al time system, given in the Gregorian date format. 

Dependency : None. 

The Al date, in modified Julian date format is the 
current independent variable for time in GMAT. There- 
fore, it is not necessary to convert the date to another 
system for this parameter. The only calculation required 
for this parameter is to use the algorithm in Sec. ?? to 
convert from Modified Julian date format to Gregorian 
date format. 


4.2.2 A IMod Julian 

Description: A IMod Julian is the epoch of an object, in 
the Al time system, given in the modified Julian date 
format. 

Dependency : None. 

The Al date, in modified Julian date format is the 
current independent variable for time in GMAT. There 
are no calculations required for this parameter. 
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4.2.3 Altitude 

Description: Altitude is the distance between a space- 
craft and a plane tangent to the surface of the body at 
the sub-satellite point. GMAT assumes the body is an 
ellipsoid. The equatorial radius, and properties of the el- 
lipsoid depend upon the particular body chosen by the 
user. 

Dependency: Central Body. 

Given: r in T\ 

Find: A 
Definitions: 

• T\ is the coordinate system in which GMAT origi- 
nally knows r 

• Tp is body fixed system of the central body selected 
by the user. 

• / is the bodies flattening coefficient 

• R is the bodies mean equatorial radius 

• 4> g d is the geodedic latitude of the spacecraft in the 
body fixed frame. 

• h is the Altitude parameter 

First we calculate <j> 9 d using the algorithm in Sec. 
4.2.21. However, to calculate h. GMAT does not convert 
to degrees, or use the modulo function. 

Then, with r expressed in Tp, we perform 


r xv - yjx 2 + y 2 

(4.93) 

e 2 = 2/ - / 2 

(4.94) 

r R 

1 a: y 1 V 

(4.95) 

CQ ^ri) y 1 — e 2 sin 2 <j> gd 


4.2.4 AOP 

Description: AOP is the argument of periapsis of a space- 
craft. The argument of periapsis is the angle between 
the eccentricity vector and a vector in the direction of 
the right ascension of the ascending node. See below for 
treatment of circular and equatorial orbits. This algo- 
rithm is adopted from Vallado. 1 

Dependency: Coordinate System. 

Given: r and v 


v - ||v|| 



Special Case: Circular Orbit 

if e < 10 -11 then, = 0.0 and return. 

Otherwise continue, 

h = rxv 


h= ||h 



Special Case: Elliptic, Equatorial Orbit 
if i < 1 0 ' 1 then, 



where e x is the first component of the eccentricity vector. 
Fix quadrant, for w: if e y < 0, then u> = 2tt — w 

Otherwise continue 

Special Case: Elliptic, Inclined Orbit 

n = [ 0 0 1 ] T x h 

_j n ■ e 

UJ COS 7] mi ' fi 

HNI 

Fix quadrant for u>: if' t z < 0, then u> = 2?r — ux 
Finally, is converted to degrees. 

4.2.5 Apoapsis 

Description: Apoapsis is the parameter used in stopping 
conditions to allow' the stopping condition algorithm to lo- 
cate the time when a spacecraft is at apoapsis. Apoapsis 
is defined as a point, along an orbital path, when the com- 
ponent of velocity, in the spacecraft position vector direc- 
tion, changes from positive to negative. The Apoapsis 
parameter is defined as the dot product of the position 
and velocity vectors. 

Dependency: Central Body. 

Given: r, v in T\ 

Find: .4 

Definitions: 

• J-\ is the coordinate system in which GMAT origi- 
nally knows r and v 


Find: ui 


r — l! r il 
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• T'i is a system with the MJ2000Eq axes, centered 
at the central body selected by the user. 

• A is the Apoapsis parameter 

if {T\ Tz) convert r and v to Tz- Then, 

A = r • v (4,97) 

4.2.6 AZI 

Description : AZI is the azimuth angle of a spacecraft, as 
shown in Fig. 4.1 using the symbol at. 

Dependency. Coordinate System. 

Given: r, v and T 

Find: aj 

AZI is calculated using the algorithm shown in Sec. 
4.1.6. There is little benefit using a routine that calculates 
only t\f and not 'tp. 

4.2.7 BdotT and BdotR 

Description: The “B” vector, B, is only defined for hy- 
perbolic orbits and is the vector from the center of mass 
of the central body, to the incoming hyperoblic asysmp- 
tote, such that the length of B is a minimum. Another 
way to say this is that B is perpendicular to the incoming 
asymptote. Let’s define S as a unit vector in the direc- 
tion of the incoming asymptote. Then, T is a unit vector 
perpendicular to S, that lies in the xy-plane of the coor- 
dinate system, Tb, chosen by the user. R is a unit vector 
perpendicular to both S and T. Finally, BdotT is the dot 
product of B and T, and BdotR is the dot product of B 
and R. The method below was adopted from work by 
Kizner. 6 

Dependency: Coordinate System. 

Given: r, v, and definition of Tb 
Find: Bn and Hr 
Definitions: 

• T\ is the coordinate system in which GMAT origi- 
nally knows r and v 

• Tb is the coordinate system in which to perform B- 
plane calculations. GMAT will place T in the xy- 
plane of Tb- Tb must have a gravitational body at 
its origin. 

• p is the gravitational parameter of the central body 
at the origin of Tb 


• B f 1 is the dot product of B and R 

• Bt is the dot product of B and T 

if the selected coordinate system does not have a ce- 
lestial body as its origin, then exit and throw an error 
message. 



Figure 4.3: Geometry of the B-Plane as Seem From a 
Viewpoint Perpendicular to the B-Plane 



Figure 4.4: The B- Vector as Seen From a Viewpoint Per- 
pendicular to Orbit Plane 

if Ti =£- Tb convert r and v from T\ to Tb 
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4.2.8 BetaAngle 

Definition: The Beta angle, j3, is defined as the angle 
between the orbit normal vector, and the vector from the 
celestial body to the sun. 

r© x v© 


e = jjejl 

„ e 

e = - 
e 

If c < 1, then the method fails and returns. 


n = 


r s © = 


Ik® x v ©! 

G© 


IU®I 


8 — sin 1 ( h • 


(4,99! 


r - ||r|| 
t; =- ||v|| 

Calculate eccentricity related information 
fv 2 — — Jr — (r • v)v 


Now let’s calculate the angular momentum and orbit 
normal vectors. 

h — r x v 

h — ||r x v|| 



• r©: Position vector of spacecraft with respect to 
celestial body, in the EarthMJ2000Eq system. 

• v©: Velocity vector of spacecraft with respect to 
celestial body, in the EarthMJ2000Eq system. 

• r.,®: Position vector from celestial body, to the sun. 


A unit vector normal to both the eccentricity vector 

and the orbit normal vector is defined as: 4.2.9 BVector Angle and BVectorMag 


h = h x e 

The following relations are only true for hyperbolic orbits: 
The semiminor axis, b, can be calculated using 

b f 

fl.yj e 2 — 1 

The incoming asymptote is defined using 



To avoid code reduplication, the magnitude and angle of 
the B vector, ||B|| and &b respectively, are calculated from 
the outputs of the B-Plane coordinates algorithm. The 
equations for j]B|| and 0 B are 

II B || = y 'Bl + Bl (4.100) 

, Br 

0 B — tan 1 ----- (4.101) 

which is implemented using atan2{Bn, Bt) 


The B-vector, B, is calculated using 
The remaining vectors, T and R are found using 



4.2.10 C3Energy 

Given: a. and jj. 

Find: C 3 


T - -- 


Sy ~S X 0 j T 

v f sfTs* 


R = S x T 

Finally, the desired quantities are found using 


C 3 — — — (4.102) 

a 

Comment: a is calculated from the satellite cartesian 
state as shown in Section 4.1.2, and p is associated with 
the specified central body. 


B t = B T 

b r - b r 


4.2.11 DEC 


if T\ -■£ Ti , convert r and v to T?, 


Description : DEC is the declination of a spacecraft, as 
shown in Fig. 4.2 using the symbol 8 . 


A = r - v 


(4.98) 


Dependency: Coordinate System. 



4.2. SIMPLE PARAMETERS 
Given: r, v and JF 
Find: 5 

Begin by converting r and v to T if necessary. Then, 
r — ||r|| (4.103) 

tf = sin _1 (^) (4.104) 

4.2.12 DECV 

Description: DECV is the declination of velocity of a space- 
craft. 

Dependency: Coordinate System. 

Given: r, v and T 
Find: 5 V 

Begin by converting r and v to T if necessary. Then, 

= s:in -1 (— ) (4.105) 

v 

4.2.13 ECC 
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Given: r, v, and coordinate system T . 

Find: ip 

Begin by converting r and v to T if necessary. Then, 

ip = cos" 1 (4.110) 

\ rv J 

4.2.15 EA 

Given: e 

Find: E 

If e > ( 1 - le -11 ) then E — 0, return. 

Otherwise, 


sm(E) = 

\/l **- e 2 sin(i/) 

1 + e cos v 
e + cos v 

(4.111) 

cos (E) = 

(4-112) 

1 + e cos v 

E = 

atan2 (sin E ; cos E) 

(4.113) 


Description: ECC is the eccentricity of an orbit and must 
be greater than or equal to zero. The eccentricity contains 
information on the shape of an orbit. If ECC is zero then 
the orbit is circular. If ECC is greater than zero, but less 
than one, the orbit is elliptic. If ECC equals one, the 
orbit is parabolic. Finally, if ECC is greater than one, 
the orbit is hyperbolic. The algorithm used in GMAT to 
calculate SMA is adopted from Vallado. 1 


4.2.16 Energy 

Description: Energy is the orbit energy. 
Dependency: Central Body. 

Given: r, v, and central body. 

Find: £ 


Dependency: Central Body. 

Given: r, v, and p. (Central Body) 
Find: e 


r = |jr || 

(4.106) 

v - || v|| 

(4.107) 

(v 2 — — )r — (r ■ v)v 


r 

(4.108) 

p 

e = j|e|| 

(4.109) 


Begin by converting r and v to a coordinate system 
with the origin equal to the central body defined by the 
user, and the MJ2000Eq axis system. Then, 


^ 2 r 


(4.114) 

(4.115) 

(4.116) 


4.2.17 HMAG 


Description: HMAG is the magnitude of the orbit angular 

4.2.14 FPA momentum. 

Dependency: Central Body. 

Description: FPA is the orbit vertical Flight Path Angle 

as as shown in Fig. 4.2 using the symbol ip. Given: r, v, and central body. 

Find: h 


Dependency: Coordinate System. 
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Begin by converting r and v to a coordinate system 
with the origin equal to the central body defined by the 
user, and the MJ2000Eq axis system. Then, 


h = l* x v 
h= llhii 


(4.117) 


(4:il8) 


4.2.21 


h ----- h 



Latitude 


(4.124) 

(4.125) 


4.2.18 HX.HY, and HZ 

Description: HX.HY, and HZ are the components of the 
orbit angular momentum vector. 

Dependency: Coordinate System. 

Given: r, v, and coordinate system T . 

Find: h x , h y , and h z 

Begin by converting r and v to T if necessary. Then, 
h — r x v — [ h x h y h z \ T (4.119) 

4.2.19 HA 

Description: HA is the orbit Hyperbolic Anomaly and is 
only defined for hyperbolic orbits. For non-hyperbolic 
orbits, HA returns a value of zero. 

Dependency: Central Body. Given: u, e 

Find: H 

If e < ( 1 + le -11 ) then H = 0, return. 

Otherwise, 


sinh (H) 

sin (z;)\/e 2 — 1 
1 + e cos v 

(4.120) 

cosh (H) 

e + cos v 
1 + e cos v 

(4.121) 

H 

! sinh H 

- *“ h 'corf, id 

(4.122) 


4.2.20 INC 

Description: INC is the inclination of an orbit in the cho- 
sen coordinate system. 

Dependency: Coordinate System. 

Given: r, v, and coordinate system mathcalF. 

Find: e 

Begin by converting r and v to mathcalF if necessary. 
Then, 

h — r x v (4.123) 


Description: Latitude is the geodetic latitude of a space- 
craft.. The geodedic latitude is defined as the the angle 
<p gc , as shown in Fig. ( ), where the sub-satellite point 
is defined by the intersceetion of a line drawn from the 
spacecraft and perpendicular to a plane tangent to the 
surface of the body. GMAT assumes the body is an ellip- 
soid, The equatorial radius, and properties of the ellipsoid 
depend upon the particular body chosen by the user. The 
algorithm in GMAT is taken from Vallado. 1 



Figure 4.5: Geocentric and Geodetic Latitude 

Dependency: Central Body. 

Given: r in T\ 

Find: <j> gc 
Definitions: 

• T\ is the coordinate system in wdiich GMAT origi- 
nally knows r 

• 'Ff is body fixed system of the central body selected 
by the user. 

• f is the bodies flattening coefficient 

• R is the bodies mean equatorial radius 

• <p g ,i is the geodedic latitude of the spacecraft in the 
body fixed frame. 
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if T\ -i- Tf, convert r from T\ to Tw. Then, 

r xy = \/x 2 + y 2 (4.126) 

Calculate the geocentric latitude to use as an initial guess 
to find the geodetic latitude 

4> g d « atan2(«,r xl ,); (4.127) 

e} — ‘If — f 2 (4.128) 

Set 5 -■ 1.0 to initialize the loop, then, 

While ( <5 > 10" 7 ) 


$gd 

(4.129) 

_ ( Re 2 sin z (b Q d 

1(4.130) 

Btan2 2 + —7==========, r xy 

\ y i - eT sin <p gd J 

! 

\4>gd - 4>’\ 

(4.131) 


EndWhile 


Dependency: Central Body. 

Given: r, ti (epoch of spacecraft in internal time system), 
and central body 

Find: O^st 

Definitions: 

• Ti equatorial inertial system (as described in Sec. 
3.4.1) of selected central body. 

• T'f is the central body’s fixed system (as described 
in Sec. 3.4.9) 

• A is the longitude of the object in Tp 

• Omst is the mean sidereal time of the central body’s 
prime meridian. 

• ti (epoch of spacecraft in internal time system) 


After convergence, 6 g d is converted to degrees, and 
converted to fall between —90° and +90° degrees. 


4.2.22 Longitude 


Description: Longitude is the longitude of an object, in 
the body fixed frame of the central body chosen by the 
user. 

Dependency: Central Body. 

Given: r, central body. 

Find: <p 

Begin by converting r to the body fixed system of the 
central body defined by the user. Then, 

4> — tan .2 1 (y, a:); (4.132) 

The calculation is completed by converting to degrees and 
setting the value to such that —180 < <j> < 180. 



Figure 4.6: Local Sidereal Time Geometry 


4.2.23 LST 

Description: LST is the local sidereal time of an object, 
with respect to the selected central body. The local side- 
real time is the sum of the longitude in the bodies fixed 
frame, and the mean sidereal time. This is illustrated in 
Fig. 4.6, where Ti is the body’s equatorial inertial sys- 
tem (as described in Sec. 3.4.1), Tf is the body’s fixed 
system (as described in Sec. 3.4.9). A is the longitude 
of the object, in this case a spacecraft, and Qmst is the 
mean sidereal time of the prime meridian. 


We begin by calculating A using the algorithm de- 
scribed in Sec. 4.2.22. The mean sidereal time Oust 
is calculated differently for Earth than for other central 
bodies. Tf the central body is Earth, then we use the 
following equations to calculate Omst- 

First, convert ti, which is the spacecraft epoch in the 
interal time system (A1 Modified Julian Date), to Tuti-, 
which is the number elapsed Julian centuries from the 
J2000 epoch. 


tuti - 21544.5 
36525 


(4.133) 
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@m ST — 67310.54841*+ 

(876600 /! ( + 8640184.812866)'i’ cm + 

0.0931042 y T1 - 6.2 x 10 _6r i$ T1 

(4.134) 


Comment: a and e are calculated from the satellite carte- 
sian state as shown in Section 4.1.2, and ft is associated 
with the specified central body. 


4.2.27 OrbitPeriod 


4.2.24 MA 

Given: v, e 
Find: M 

If e < ( 1 - le~ n ) then calculate E using algorithm in 
Sec. 4.2.15. Then M is calculated using 

M = E — e sin E (4.135) 

Note: E must be expressed in radians in the above equa- 
tion, and results in M in radians. 

If e > ( 1 + le~ n ) then calculate H using algorithm in 
Sec. 4.2.19. Then M is calculated using 


Given: a, and u 
Find: I' 

If a < 0, then 1 ' = 0, return. 

Otherwise, 

la 3 

T — 2ttW — (4.140) 

V P 

Comment: a is calculated from the satellite cartesian 
.state as shown in Section 4.1.2, and u is associated with 
the specified central body. 


M = eCmh H - H (4.136) 4.2.28 PercentShadow 


Note: H must be expressed in radians in the above equa- 
tion, and results in M in radians. GMAT outputs MA in 
degrees. 

If neither of the above conditions are satisfied, M ----- 0, 
and output “Warning: Orbit is near parabolic in mean 
anomaly calculation. Setting MA = 0” . 


4.2.25 MHA 

4.2.26 MM 

Given: a, e, and ft 
Find: n 

The orbit is considered either circular or elliptic ( both 
orbit types use the same equation to calculate n) if e < 
1 — le~ lL . In this case the mean motion, n, is calculated 
using 

» + (tm) 

The orbit is considered hyperbolic if e > 1 + le ll . In 
this case the mean motion, n, is calculated using 

n=W-4 14.138) 

V a d 

If neither of the above two conditions are met, the mean 
motion is calculated using 


The PercentShadow parameter calculates the percentage 
of the apparent solar disk that is in view from the per- 
spective of a spacecraft. The algorithm used in GMAT 
was adapted from Montenbruck 7 pgs. 80-83. 



Figure 4.7: Shadow Geometry 


• Rq = Radius of the Sun 

• Rb -= Radius of occulting body 


RL 


Apparent radius of the Sun 


n — 2 v /q 


(4.139) 
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• R ! g ----- Apparent radius of occulting body 

• r v, — Vector from central body to Sun 

• tb — Vector from central body to occulting body 

• r = Vector from central body to s/c 


We begin by calculating the vector from the occulting 
body to the spacecraft, s, using 

s = r - rfj (4.141) 

and the vector from the occulting body to the sun, s®, 
using 

s Q = r® — t b (4.142) 

(Note that when the occulting body is the same as the 
central body, s = r, and s® = r®) 

Next we calculate the apparent radius of the Sun and 
occulting body' using 




Figure 4.8: Occultation Geometry in Calculation of Per- 
centShadow 



We can calculate the apparent separation of the two bod- 
ies, V , using 


D ! 



( ~s r (r® -r) A 

\ s ll r © — r ll J 


(4.145) 


If O' > Rq + R'f S , then the spacecraft is not in the 
body’s shadow and 


p = 0: 


(4.146) 


4.2.29 RA 

Description: RA is the right ascension of a spacecraft, as 
shown in Fig. 4.2 using the symbol A. 

Dependency: Coordinate System. 

Given: r, v and T 


If O' < R' b — RL, then the spacecraft is in full shadow 
and 

p = 100; (4.147) 


If neither of the above conditions are met, the space- 
craft is in partial shadow. 

li\R' a -R' B \ < D' < R' a +R! B , then we can calculate the 
percentage of shadow by calculating the area of overlap, 
A, of the two apparent disks as shown in Fig. 4.8. 



and 

C3 = y/tfS~c ? (4.150) 


Find: A 

Begin by converting r and vtof if necessary. Then, 
A = tanJ 1 (y,x) (4.153) 

4.2.30 RAV 

Description: RAV is the right ascension of velocity of a 
spacecraft, as shown in Fig. 4.2 using the symbol A„. 

Dependency: Coordinate System. 

Given: r, v and T 

Find: A„ 

Begin by converting r and v to T if necessary. Then, 
A„ = tan 2 ' 1 (t>, / ,r; a; ) (4.154) 
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4.2.31 RAAN 


4.2.34 RMAG 


Description: RAAN is the right ascension of the ascending 
node as shown in Fig. 4.1 using the symbol Q. 

Dependency : Coordinate System. 

Given: r, v, and coordinate system T . 

Find: e 


Description: RMAG is the magnitude of the spacecraft’s 
position vector. 

Dependency : Central Body. 

Given: r and central body. 

Find: r 


Begin by converting r and v to F if necessary Then, 
h — r x v (4.155) 


Begin by converting r to a coordinate system with the 
origin equal to the central body defined by the user, and 
the MJ2000Eq axis system. Then, 


h :=== ||hjj 

(4.156) 

n = | 0 0 1 j T x h 

(4.157) 

3 

1! 

jr 

(4.158) 

*— “(t) 

(4.159) 

if (i > IQ" 11 ), then 


Q = cos~ x ( — 

\ Tl / 

(4,160) 


Fix quadrant for C: if n y < 0, then Q = 27r — 0 
if (i < 10 —li ), then 

Q :: 0 (4.161) 


r — |jr|| (4.164) 

4.2.35 SemilatusRectum 

Description: SemilatusRectum is the orbit semilatus rec- 
tum, which is the magnitude of the position vector, when 
at true anomaly of 90°. 

Dependency: Central Body. 

Given: r, v, and ji. (central body). 

Find: p 

Begin by converting r and v to a coordinate system 
with the origin equal to the central body defined by the 
user, and the MJ2000Eq axis system. Then, 


4.2.32 RadApo 

Given: a, and e 
Find: r a 

if 1 — e < 10 -12 then r a = 0. Note, this means that for 
parabolica, and hyperbolic orbits, GMAT outputs a value 
of zero for RadApo. Otherwise, 


h = r x v 

(4.165) 

h = ||h|j 

(4.166) 

h 2 

P = — 

(4.167) 

T 


4.2.36 SMA 


r a = o(l + e) (4.162) 

Comment: a and e are calculated from the satellite carte- 
sian state as shown in Section 4.1.2. 

4.2.33 RadPer 

Given: a, and e 
Find: r p 


Description: SKA is the semimajor axis of an orbit. The 
SMA contains information on the size and type of an 
orbit. If the SMA is positive, the orbit is elliptic. If 
the SMA is negative the orbit is hyperbolic. The SMA 
is undefined for parabolic orbits. The algorithm used in 
GMAT to calculate SMA is adopted from Vallado. 1 

Dependency: Central Body 

Given: r, v, and p, (Central Body) 

Find: a 


r p ----- a(l — e) (4.163) 

Comment: a and e are calculated from the satellite carte- 
sian state as shown in Section 4.1,2. 


r — Ill'll 

* = IMI 



(4.168) 

(4.169) 

(4.170) 
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1 


V V K) 1. ,l\ 


o. 


s_. 


if 11 -- el > 10 30 . then 


JL 

2f 


(4.171) 


otherwise, report error and return. Error: “Warning: Or- 
bit is near parabolic and SMA is undefined”. 


4.2.37 TA 

Description: TA is the orbit true anomaly as shown in 
Fig. 4.1 using the symbol u. 

Dependency: Central Body. 

Given: r, v, and coordinate system T . 

Find: u 

Begin by converting r and v to T if necessary. Then, 


h — r x v 

h=\M 

n — [ 0 0 1 ] T x h 
n — ||n|| 

r — l! r ll 
« = HI 

(v 2 ■••• — )r - (r • v)v 
r 

e - li e il 


l — cos 


4.2.38 TAIMod Julian 

Description: TAIModJulian is the epoch in the TAI time 
system, expressed in the modified Julian date format. See 
Sec. 2.1.1 and 2.2.1 for more details. 

Dependency: None. 

Given: A1 (epoch in the internal, Al time system). 

Find: TAI 


To convert from A1 to TAI we use the following equa- 


tion 


TAI = .41 - 0.0343817sec 


(4.184) 


4.2.39 TTModJulian 


(4.172) 

(4.173) 

(4.174) 

(4.175) 

(4.176) 

(4.177) 

(4.178) 

(4.179) 

(4.180) 


Description: TTModJulian is the epoch in the TT time 
system, expressed in the modified Julian date format. See 
Sec. 2.1.3 and 2.2.1 for more details. 

Dependency: None. 

Given: A1 (epoch in the internal, A1 time system). 

Find: TT 


To convert from A1 to TT we use the following equa- 


tion 


TT = .41 - 0.0343817sec + 32.184sec 


(4.185) 


There are three special cases, and they are treated differ- 
ently. 

Special Case 1: Elliptic Orbit 
if (e > ID" 11 ), then 

y = cos -1 f — — (4.181) 
\ er I 

Fix quadrant for v: if r • v < 0, then v =- 2tv — v 

Special Case 2: Circular, Inclined Orbit 

if (e < 10" n ) and ( i > 10 — 1 1 ), then 

v — cos -i ( 'j (4.182) 

Fix quadrant for u: if r z < 0, then v = 2rr — v 

Special Case 3: Circular, Equatorial Orbit 

if (e < 10~ n ) and (i < 10"" 11 ), then 

v = cos" 1 (— ) (4.183) 

\ r / 

Fix quadrant for u: if r v < 0, then v —2r — v 


4.2.40 TTGregorian 

Description: TTGregorian is the; epoch in the TT time 
system, expressed in the Gregorian date format. See Sec. 
2.1.3 and 2.2.2 for more details. 

Dependency: None. 

Given: Al (epoch in the internal, A1 time system). 

Find: TT 

To convert from Al to TT we use Eq. (4.185). Then, 
knowing the epoch in the TT time system in the modified 
Julian date format, we use the algorithm in Sec. 2.2.1 to 
obtain the Gregorian date. 

4.2.41 Umbra and Penumbra 

The Umbra and Penumbra parameters are used to de- 
termine if a spacecraft is in the shadow of a celestial 
body. The algorithm used in GMAT is adapted from 
Montenbruck ' pgs. 80-81. For both functions, if the value 
is less than 1, then the body is in shadow, if the function 
is greater than 1, then the body is not in shadow. 
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Figure 4.9: Geometry of Umbra, and Penumbra Regions 


For definitions of see Sec. 4.2.28. 


£ = - aTs © 
s © 

(4.186) 

Si-4 

1! 

1 

'nj 

(4.187) 

Rq + Jig 
sin a p — 

(4.188) 

Rq — Rg 

esmOtt — 

S-0 

(4.189) 

The radii of the umbra and penumbra cones, 
at distance t, are respectively 

t p and r u , 


4.2.42 UTCMod Julian 

Description : UTCMod Julian is the epoch in the UTC time 
system, expressed in the modified Julian date format. See 
Sec. 2.1.2 and 2.2.1 for more details. 

Dependency : None. 

Given: Al (epoch in the internal, Al time system) . 

Find: UTC 

To convert from Al to UTC we use the following equa- 
tion 

UTC = A1- 0.034381 7sec - A AT (4.196) 


r v = tan a p [1 + --- 

\ sin a p J 

(4.190) 

r v = tana,, ( £ — — -2— j 
\ sin a u J 

(4.191) 

Finally, if £ > 0 


[ 

II 

(4.192) 

d u — d \'^u\ 

(4.193) 

If £ > 0 d u < 0 and r u < 0, then the object is 
umbral eclipse region. 

in the total 


If £ > 0 d u < 0 and r u > 0, then the object is in the 
annular umbral eclipse region. 

If £ < 0, then the object is on the day side of the occulting 
body and is not in shadow and 


The default is to read A AT from the file named tai- 
utc.dat. A AT is the accumulated leap seconds since Jan. 
1961. 

4.2.43 VelApoapsis 

Given: a, e, and /< 

Find: v a 

If e > ( 1 - le -12 ) then v a = 0. 

Otherwise. 



d p — j d — r p \ 
d,. = | d- \r u \\ 


(4.194) 

(4.195) 


Comment: a and e are calculated from the satellite carte- 
sian state as shown in Section 4.1.2, and ji is associated 
with the specified central body. 
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4.2.44 VelPeriapsis 

Given: a, e, and ji 
Find: v p 


Iterate On: E n+ i = E n + 
Until: \E n +i — E n | < le~ 8 


M — E n + e sin E n 
1 — e cos E n 


Finally we convert the eccentric anomaly to the true anomaly 
using the algorithm given in sec. 4.3.2 



Comment: a and e are calculated from the satellite carte- 
sian state as shown in Section 4.1.2, and jx is associated 
with the specified central body. 


Hyperbolic Orbit Case 

If e > 1 then use the following algorithm: 

We begin by choosing the initial guess for the hyper- 
bolic anomaly. The initial guess depends on the value of 
the mean anomaly and the eccentricity: 

If e < 1.6 


4.2.45 VMAG 

Description: VMAG is the magnitude of the spacecraft’s ve- 
locity vector, when the velocity is expressed in the chosen 
coordinate system. 

Dependency: Coordinate System. 

Given: v and coordinate system T. 

Find: v 

Begin by converting v to coordinate system T if nec- 
essary. Then, 


If ( - 7r < M < 0 ) or M > 7T 
11 — M — e 

Else 

1-1 = M + e 

End 


Else 

If (e<3.6& \M \ >n) 
H = M — sign (Af)e 

Else 

fl-fu 

End 

End 


V ■ jjvjj = yjv% +V% + v'i 


4.3 Other Calculations 


(4.199) Iterate to determine the Hyperbolic Anomaly: 

M + 11 n — e sinh H n 


Iterate On: H n+ i = H. n + 
Until: | H n +i — H n \ < le'" 8 


e cosh H n — 1 


4.3.1 MA to TA 


Convert the hyperolic anomaly to the true anomaly using 
the algorithm given in sec. 4.3.3 


Description: This algorithm shows how to calculate v 
given M and e and is taken from Vallado. 1 

Given: M. and e. 

Find: v 

The algorithm is different for elliptic and hyperbolic 
orbits. Let’s first look at what happens for elliptic orbits. 

Elliptic Orbit Case 

If e <= 1 then use the following algorithm: 

Determine initial guess for the Eccentric anomaly 


If (- 71 

- < M < 1 

E 

= M - e 

Else 


E 

=■■ M + e 

End 



Iterate to determine the eccentric anomaly: 


4.3.2 EA to TA 

Description: This algorithm shows how to calculate v 
given E and e and is taken from Vallado. 1 

Given: E and e. 

Find: v 


y/\ — (E sin(iv) 
Sill v — — — 

(4.200) 

1 — e cos b 

cos E — e 

(4.201) 

cos v ~ — 

1 — e cos h 


i / = atan2(sin u. t cos u) 

(4.202) 



vvryr a 
V V KJ i l.\. 


4.3.3 HA to TA 


Description : This algorithm shows how to calculate v 
given H and e and is taken from Vallado. 1 

Given: 11 and e. 


CHAPTER 4. CALCULATION OBJECTS 

This coordinates system is illustrated in Fig. 4.11. 
The locations of the libration points in the rotating coor- 
dinate system can be found by calculating the values of 
7 that solve the following equations: 


if - (3 - /A) 7 f + (3 - 2 /t *) 7? - A 7? 

+ 2/i*7i - /A = 0 (For LI 


^ (4.206) 


(4.203) 


(4.204) 

(4.205) 


\/e 2 — 1 sirili ( // ) 

smi/ ---- ---- 

1 — e cosh H 

cosh H — e 
008 " = 1 - e cosh H 

v = atan2(sin u, cos n) 


4.4 Libration Points 


We begin by assuming that the planets move in circu- 
lar orbits about the sun, and the mass of a spacecraft is 
negligible compared to the mass of the planets. For il- 
lustrative purposes, lets consider the Earth and its orbit 
about the Sun. In this case, the libration points are lo- 
cations in space where a spacecraft, will stay fixed with 
respect, to the Earth and Sun. Figure 4.10 shows a simple 
illustration. We see the Sun, the Earth’s position with 
respect to the Sun, and the Libration points Lj and Lz 
at two different epochs. Notice that at. f 3 , the points L r 
and L '2 are on the Earth-Sun line. At a later time, £2, 
although the Earth has moved with respect to the sun, 
L] and L-z still lie on the Earth-Sun fine. 

The preceding example gives a brief qualitative de- 
scription of two of the Earth-Sun libration points. In 
general, there are five libration points for a given three 
body system. To determine the locations of the libration 
points, it is convenient to work in a rotating coordinate 
system rather than the inertial system shown in Fig. 4.10. 
The system we use is constructed as follows: 

• Define the primary as the heavier of the two bodies, 
the secondary as the lighter. 

• Define the coordinate system x-axis as the axis point- 
ing from the primary to the secondary. 

• Define the y-axis to be orthogonal to the x-axis in 
the plane of the secondary’s motion about the pri- 
mary, pointing in the direction the secondary moves 
about the primary. 

• Define the z-axis orthogonal to the x and y axes to 
form a right-handed system. 


7! + (3 - in 7! + (3 - 2 p*) ji - nil 


- 2/A 72 - 


= 0 (For L2) 


7f + (2 + /A) 7 3 + (1 + V) 7;s — ( 1 — n) 7 f 
- 2(1- /A) 73 -• (1 - /A) = 0 (For L3) 


where 


m-i + m 2 


(4.207) 


(4.208) 


(4.209) 


Equations (4.206)-(4.208) do not have exact analytic 
solutions.Szebehely 8 notes that they are most easily solved 
using an iterative method with the following as the initial 
guesses: 


71 = 72 = 


\3 (1 - /A) 

/3 — 1 


(4.210) 

(4.211) 


GMAT uses the Newton-Raphson method to solve for 
the roots of the equations by iterating on 


7(1 + 1) = 7 (i) 


F(y(i)) 
A'(7 W) 


(4.212) 


until the the difference (7 (i + 1) — 7 (i)\ < 10 8 . The 
derivative F'( 7) for each libration point is shown belov/. 


P'{ 7) = oof - 4 (3 - /T)7 ? + 3 (3 - 2/A) 7? 


- 2/t*7i + 2/A (Fbr LI) 


(4.213) 

I A ' 


^(7) = 5 7 | + 4 (3 - /A ) 7 | + 3 (3 - 2/A) 7 2 

- 2/4*72 - 2/4* (For L2) 

^( 7 ) = 073 + 4 (2 + /A) 7f + 3 (1 + 2/A) 7 f 

- 2 (1 - /A) 73 ~ 2 (1 - /A) (For L3) j 

We now need to redimensionalize the results found 
in the rotating system, and perform the necessary trans- 
formations to obtain the results in the M.J2000 system. 
Let’s assume that r s , v.,, and a., are the position, veloc- 
ity, and acceleration vectors respectively of the secondary 
body, with respect to the primary body, expressed in the 
FK5 system. Then, the position of the i th libration point 
can be expressed in the rotating system with the origin 
centered on the primary body as 


r s i J'i Vi 


• Place the origin at. center-of-mass of the system. 


(4.216) 
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Figure 4.10: Geometry of Libration Points 


Table 4.7: Location of Libration Points in RLP Frame, 
with the Origin at the Primary Body 


Point 

a>Position 

’./-Position 

LI 

1 — 7i 

0 

L2 

1 + 72 

0 

L3 

-7s 

0 

L4 

1/2 

\/3/2 

L5 

1/2 

-\/3/2 


where 

r. = ||r.|| (4.217) 

The velocity of the i th libration point can be expressed 
in the rotating system with the origin centered on the 
primary body as 

v* - \xi yi Of (4.218) 

r s 


and 


where 


and 


f ~ 

X: 

yi 

*A 

X2 

in 

h 

x-i 

y-3 

h J 


r s 



r a 


. r s x v., 

z = 7 t 7 

jr s x v s i 


y = z x it 


Vs 

r s 


l's 

— r s • v,) 


r s x a s 
!|r 5 XV, || 


Z 

||r s x v s jj 


(r, x a s • z) 


(4.220) 


(4.221) 


(4.222) 


(4.223) 


(4.224) 


(4.225) 


Now we have the redimensionalized position and velocity 
vectors of the libration point in the rotating coordinate 
system defined by the motion of the secondary body with 
respect to the primary body. To determine the position 
and velocity vectors in the FK5 system, with the origin 
located at the primary body, we need to determine the 
rotation matrix and its derivative as follows: 


R/ 




y 3 h) 


(4.219) 


y = zxx + zxx (4.226) 

GMAT currently assumes that the terms r, x a* are zero. 

We finally arrive at the position of the Libration Point 
in the FK5 system with the origin at the primary by per- 
forming the calculations: 


r = R/V 
v = R/V + R'V 


(4.227) 

(4.228) 
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Figure 4.11: Location of Libration Points 


4.5 B ary center 

The barycenter of a system of point masses, Tb, is also 
called the center of mass. If we have a system of n bodies, 
and we know the position of the i th body with respect to 
a common reference system, then we can calculate the 
barycenter of the system using 

11 

y_ . niiTi 

tb - — (4-229) 

y mi 

i — 1 


Similarly, we ca.n calculate the velocity of the barycenter 
using the following equation 


n 

y mvi 

i—1 


n 



i-1 


Vfe = 


(4.230) 



Chapter 5 

Dynamics Modelling 


One of the fundamental capabilities of GMAT is to 
model the motion of spacecraft in many different flight 
regimes. The flight regimes, such as low Earth, or Li- 
bration Points, are determined by which forces and per- 
turbations dominate the dynamics. In this chapter we 
present how GMAT models the dynamics of spacecraft in 
motion. We discuss how GMAT calculates many different 
types of forces including mnlitiple non-spherical gravity- 
perturbations. third-body effects, and atmospheric drag 
among others. We being by looking at the general form 
of the equations of motion. 


d 2 r 

di 2 


”3 r+Wfc+Ggm* 

fc-I 

Mi 




+£(v#. 


k= 1 

kj=i 


_ \ rh.. dr 

V^ +— ^7 
J m at 


1 2 

- 2^ -^7 


V re / “1“ 


I^srCrAq 


m. 


-r s0 


(5.3) 


Description 

Central Body Point Mass 


Term 



5.1 Equations of Motion 


Central Body Direct Nonspheri- 
cal 


5.1.1 Orbit State Equations 


Direct. Third Body Point Mass 



Let’s begin by defining the position and velocity of a 
spacecraft with respect to the central body' of integration 
as rand v. From Newton’s Second Law we know that 


Indirect Third Body Point Mass 


nt 






(5.1) 


Third Body Direct Nonspherical 


£(v«.) 


which says that the mass, times the acceleration, is equal 
to the sum of the forces. Solving for the acceleration gives 
us the second order differential equation 


d 2 r 
t iH 



(5.2) 


Third Body Indirect Nonspheri- 
cal 



Spacecraft Thrust 


771,5 dr 

m dt 


Atmospheric Drag 


1 2 C d A 

nP^rel __ X rel 


CMAT has the capability to model many different types 
of accelerations experienced hv spacecraft in orbit. If we 
include all of the possible forces GMAT can model in the 
summation on the left hand side of Eq. (5.2). then we 
would have 


Solar Radiation Pressure 


^srCrAq „ 

r sG 


In general, Eq. (5.3) does not have an analytic solution 


■57 
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so CM AT uses numerical integration to find approximate, 
although very accurate, solutions. GMAT uses first order 
numerical integrators, so we must take the three second 
order differential equations in Eq. (5.3) and convert them 
to six first order equations. So, we define a new variable 
x such that 


5.1.3 Multiple Spacecraft Propgation and 
Coupled Propagation of the Equa- 
tions of Motion 

5.2 Force Modelling 


x = [r r v J ] T — [x y z x ij 

then taking the derivative we arrive at 
x— [r r v / ] T — [i ij z x ij 


z} 1 (5.4) 

if (5.5) 


5.1.2 State Transition Matrix Equations 


where 




A &{t,t 0 ) 



(5.6) 

(5.7) 


5.2.1 Non-Spherical Gravity 

GMAT integrates all spacecraft equations of motion us- 
ing the Earth’s Mean J2000 axis system. However, the 
user can choose central bodies other than the Earth as 
the origin of the coordinate system of integration. Grav- 
itational forces are conservative and only a function of 
position. To calculate the gravitational force due to a 
non-spherical body, we need to determine the position of 
the spacecraft in the body fixed frame Tf- However, the 
equations of motion are expressed in terms of the position 
of the spacecraft in the inertial frame. 

We know from dynamics that the acceleration in an 
inertial frame can be calculated using 


subject to the initial conditions 


a c6 - V<7 


(5.16) 


*(t 0 ,t 0 ) 1-6x6 


(5.8) 


If we define x as 


x — 


r 

v 


then. 



Now we can write A as 



( dv 

dv \ 

dx 

dr 

dv 

dx ~ 

da 

da 


\ 9r 

Ihs / 


(5-9) 


(5.i0) 


(5.11) 


For convenience, lets use the following notation 


A = 

dv 

(5.12) 


dr 

B = 

dv 

(5-13) 

C = 

da 

(5.14) 

D = 

dr 

da 

(5.15) 

dv 


where U is the gravitational potential. The potential for 
a nonspherical body comes from the solution to Laplace’s 
equation: 

V Z U - 0 (5.17) 

The solution to this equation is most easily' expressed in 
spherical, body-fixed coordinates because it allows for a 
convenient separation of variables. 

In spherical coordinates the gradient of the gravita- 
tional potential is 


L dU 1 dU 

- H - U A 

' o<p r COS <p OA 


(5.18) 


We see that there axe two singularities in Eq. (5,18). The 
first is when r = 0, which is a nonphysical case and we 
will not discuss it further. The second singularity occurs 
when <j> -■=■ ±90°. Pines 9 developed a uniform expression 
of the gravitational potential that avoids the singularity 
a,t the poles: 


U 


00 f tL \ n n 

1 + /L ( I 'll Anrn ip)\Gn m cos(m\) cos’" $ 
n— 1 ' y> in — 0 


+ S nm sin(rn\) cos’" <f>] 


(5.19) 


Examining this form of the potential it is easy that there 
is not a. singularity at the poles when taking the gradient 
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in spherical coordinates. Pines rewrites Eq. (5.19) as 


l/=£ 


Ro 


i+E 

n— 1 


nrn (v., 1 L f) 


(5.20) 


where CU and S' nm are the gravitational coefficients, s, 
t, and « are given by 

s = x/r, t ~= y/r, u ~ zjr — sin <?■ 

and r m (s, t ) and i, n (s, t) are calculated using the recursive 
relationships 

r 0 = 1, n = s, i 0 ==■ 0, U = t 
f m ~ Sf m~\ tlm —\ , i m = 1 T tfm—1 

The coefficients j4 ntn («) are called “derived” Legendre 
functions and are given by 


The derived Legendre polynomials are normalized using 

A nm — N nm A n7n (5.27, 

where A rem are the normalized Legendre polynomials. 
Lundberg 10 showed that there are several recursive al- 
gorithms to compute A n - m but that only two are stable. 
GMAT uses the following algorithm to recursively calcu- 
late the derived Legendre polynomials 


■1?: y/i 


(2n + l)(2n • 


1/2 


An- 


l,?n 


(n — m)(n + m) 
f (2n + l)(n — m — l)(n + m — 1) 

L 


(2n — 3) (n + rn) (n — m) 


1/2 


The recursive algorithm is started using 
Au = V3 cos <t> 

./Li.yy — COS 


d m 

A nm (U) — 

where we know from Rodrigues’ 10 formula that 

1 d 


l‘2n +l j 

/ - An — ■ 1 ,n~ 1 


2 n 


A n ~2.m 

(5.28) 

(5.29) 

(5.30) 


(5.21) 


and 


a»W = a W = - ^ (5.22} 


«,„(«) - (1 - C,(») (5.23) 


The above equations are normalized using Eq. (5.27) and 
used in 

The acceleration due to nonspherical gravity can be 
written as 


a ff 


For numerical reasons it is useful to normalize some 
of the terms in the potential function, U. By normal- 
izing the spherical coefficients and the derived Legendre 
polynomials we can improve the stability of recursive al- 
gorithms used to calculate the Legendre polynomials and 
improved numerical problems. We use the nondimension- 
alization approach and described by Lundberg. 10 Lund- 
berg chooses the normalization factor so that the normal- 
ized spherical harmonics C nm and S nm will have a mean 
square value of one on the unit sphere. The normalized 
Legendre functions, P nm , are defined so that the product 
of the spherical harmonic coefficients and the correspond- 
ing Legendre functions remain constant, or 


dU s dV t dU udU 
dr r ds r dt r du 
( 1 3V 1 dU 1 dU\ r 

lr 3s r 8t r du ) 


(5.31) 


To simplify the partial derivatives in Eq. (5.31), Pines 
defines some intermediate variables as follows 


Po = p/r 

Pi = PPo 

Pn = ppn— i for n > 1 


(5.32) 


FnrrtCn 


I iLinCnir 


Pirn ‘dura — I nm^nm 

(5.24) 

GMAT uses the normalization factor N nm given by 


An 


(n — rn)!(2n + 1)1 1 1/2 


(n + m ) ! 


(5.25) 


The non-dimensional spherical harmonic coefficients and 
Legendre functions are 

P = N P C — ~ n,rb § nm 

1 in n — 1 y nm l mn nm — *»'n 


Using Limdberg’s nondixnensionalization approach, we can 
write 

Uyym(A, f) — G nm r., n (s , t) -{- S nm i m ( s •. 

~ Cnm'fm -- 1 (*b f) T 1(*5 I) 

Ani nf^yf) “ R-ninf m — 1 ( s > ^) Pnm^ J m— 1(^> 0 
C n, rl (s , t) — C nm r m —2(.S, t) + Smnim— 2( s j I ) 
ILvm'Oi t) — S n! nf m—2{$ . 1) Cnm^m— 2(^,1) 

(5.33) 

The partial derivatives in Eq. (5.31) can be written as 

t 3U u dU 
r dt r du 


dU 

dr 


s_dU_ 

r ds 

n 


A„ 


a; 


(5.26) 


" u ' y ^ 1 , m + 1 -^7/.+ 1, m+ 1 ^mn 

n—0 ^ m— 0 


(5.34) 
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1 fiTT oo 'W 

-■g- = E ~n~ E A n , n (u)mE n 

r ds , tr , ^ “ 0 


1 <91/ _ r- v ,0 7l+ l <r~. v . r, 

r m ~ a 2 E A nm ( u ) mK , 


n — 0 


m=0 


where 


, where C 9 is a symmetric matrix with components given 


(5-36) 


- r\TT oc n, 

- E = E ^ E ^ m +iA n , m +i(u)D nm (5.37) 

n=Q ^ m— 0 


Cll 

<42 

<4.3 

C 22 

c 23 

C 33 

Note that 


tin 4* 5 ^tl 44 -{- 0 - 4 yV ~f- 2 5 (Tl 14 (5.48) 

C 21 = <2-12 *T* 6‘ia 44 + 51224 4~ *a 14 (5.49) 

C 31 = Cf 13 4~ £ 144 X 44 + 5034 4“ tiOi 4 (5.50) 
-On + < 2 fl44 + 04 /r + 2t«24 (5.51) 

C 32 = <223 ~r tu<Z 4 4 + 1-IC124 4" £#34 (5.52) 

033 4- w"a 44 4- a 4 /r -+- 2 * u * 0,34 (5.53) 


c«,m+i — [(n - m)(n 4- m + 1)] 1/2 

(n + m 4- 2)(n + rn + 1) "1 1/2 


<41+1,771+1 


a 4 


317 sdU tdU udU 


(2ro -f 3)(2n + 2) 

To calculate the nonzero portion of the sensitivity ma- 
trix, we begin by calcluting the following 9 terms: 


dr r 9s r dt r du 
and is given in Eq. (5.34). 

5.2.2 /.'-Body Point Mass Gravity 


(5.54) 


Oil 


«12 = 


Ol3 


«14 


E 

Pn+2 

Rl 

^ m(m - 

-114 G 

x J- l *nm nm 

(5.38) 
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(5.39) 
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^ ^ ^<^?i,m+ 1 -^i,m+l &nm 
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(5.40) 
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y y \ /- '-'■ V ,. . - . 7 — } i - t 1 1 , 7 -{ i (^.4 - ^ 


71 .-O m— 0 


«23 = 5 Z 53 mc ^"*+lAM»+l2Wm (5.42) 

n=0 ® m=0 
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^ ^ ^ ^ 4 4 >.4 3 j 

71= 0 ® 771=0 
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OO 71 
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^ y / ,, <4l, 7X1+2 -^71, 771+2-^ 71777 , (5.44 y 

77=0 ^ 771=0 
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^ ] D 2 y v <41+1,771+2-^71+1,771 + 2^77,7^5.40; 

'71=0 771=0 

1244 = 

OO 71 

^ y V <^n+2, 771+2-^71+2,771+2^71771 (5.46y 

n =0 ^ m =0 

where 



Cn+l.m+2 = Cn+l.m+1 [(« - m) (« + TO + 3)] 1/2 


c It ,m+i [(« — m - l)(n -i- m + 2’ 


, 1/2 
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(n T m + 4) (n + m + 3) 


(2 n + 5)(2n + 4) 


L/2 


n _ da a 

!> “ Q t 


(5.47) 


The gravitational perturbation due to n point masses is 
well know. However, we will derive the governing differ- 
ential equation here, as well as the componenents of the 
sensitivity matrix. Let's begin by defining some notation 
referring to Fig.5.1. Assume the body is the central 
body of the integration. 



Figure 5.1: N-Body Illustration 


• r« is the position of the spacecraft with respect a 
hypothesized inertial frame. 

• Tj is the position of the central body with respect 
a hypothesized inertial frame. 

• if. is the position of the k th gravitational body with 
respect a hypothesized inertial frame. 


Finally, 
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• r is the position of the spacecraft with respect to 
the central body of integration (j ih body). 

• rfc is the position of the k th gravitational body with 
respect to the central body. 

We need the governing differential equation that de- 
scribes the motion of the spacecraft with respect to the 
central body. However, we know that we must apply New- 
ton’s 2nd Law in an inertial frame. So, we begin by defin- 
ing the relative position of the spacecraft with respect to 
the central body. From inspection of Fig. 5.1 we see that 

ij + r — r 5 (5.55) 

By reordering and taking the second derivative with re- 
spect to time we obtain 

r = r s — ij (5.56) 

We can apply Newton’s 2nd Law to the spacecraft and 
obtain 


We can break down the acceleration in the equation above 
into three physical categories. The first term is the ac- 
celeration on the spacecraft due to a point mass central 
body. The second type of terms are called direct terms. 
They' account tor the force of the k th body on the space- 
craft. The third type of terms are called indirect. They 
account for the force of the k th body on the central body. 

Let’s look at the contributions to the sensitivity ma- 
trix due to point mass perturbations. We notice that a pm 
is not a function of velocity. So. 

Npj.. — Dpm — 0;)x3 (5.63) 

We also know that 

Bpm Isx3 (5.64) 

This leaves C pm as the only noil-trivial term for point 
mass gravitational effects. Let’s look first at the deriva- 
tives of the point mass term. We can use the vector iden- 
tity in Eq. (12.4) to arrive at 


= ’^TF k = G^T 


k—l 
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m s rrik 
In “ rl 


■ (rjh — r) (5-57) 




I 3 3 



^•5 


15.65) 


where (r k — r) is a vector from the spacecraft to the k th 
body, m s is the mass of the spacecraft, and m k is the 
mass of the k th body. We can write r s as simply 


Similarly, applying Eq. (12.4) to the direct terms we 

see that 


nik 


's = (?Et ~r) 

“1 IF* - r !l 


(5.58) 
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dr 


k= 1 


\k+j 


n-r 


We can apply Newton’s 2nd Law to the j th body and 
obtain 


Grn s mj v'' m i m k 
. — i r + C > -—s-r fc 


5 N- r ll 3l3 + 3 S^ 

kjij 
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- r) Ip 


k—l 

kri 


F* 


(5.59) 


(5.66) 


where the first term is the influence of the spacecraft on 
the central body, and the second term is the influence of 
the k point, mass gravitational bodies. We can write r,- 
as simply 

71 

mk 

Ij 


Finally, the derivative of the indirect terms are zero and 
we have 

flj rr T 

t * jrl ! K — y I 3 T 5 i‘ .j — 


Gm s ^ \ - 


k—l 

kn 
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(5.60) 


ft* T , / (rfc ~r)(r fc -r) J 


-ElTli3 l3 + 3 E 
■ - ll r * — r lh ^ 


Substituting Eq. (5.58) and (5.60) into (5.56) we get 
r =GT mk 


k—l 

kj=j 


Mfc 


k—l 

kf-j 


(||r fc -r) |p 


Z—J 

k—l 


Ffc - r 


Gm. 

r a 


i ( r * - r ) Or 1 - G E i^ni3 r * 


Finally, collecting terms yields 


mk 

f=r imp 

kr-j 

(5.61) 


(5.67) 


Combining similar terms we can express the result as 
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5.2.3 Atmospheric Drag 

1 ■) G dA „ . . 

= -7iP°rer v rt: i (iJ.69) 

l TO,, 

where 

v re i = v- w 8 xi (5.70) 

where Ug is the Earth’s angular velocity vector in the 
FK5 system. 

GMAT does not currently support calculating the STM 
using drag. The components of the sensitivity matrix A 
contain derivatives of the atmospheric density with re- 
spect to position. These derivatives are non trivial for 
most density models and are not currently included in 
GMAT. 


• The user can provide a value for use in the solution 
of the equations of motion. 

The body should store the users original input for the 
state and epoch, and the state and epoch calculated at 
the last request for state information. Then, when the 
next request is made for state information, the epoch and 
state from the last request are used as the input state for 
next calculation. 

5.3.3 Atmospheric Density 

28 K p + 0.03e Kp = A p + 100 (l - e (- n - 08 A>)^ (5.77) 


5.3.4 Jacchia Roberts 

MSISE-90 

A. E. Hedin, Extension of the MS1S Thermospheric Model 
into the Middle and Lower Atmosphere, J. Geophvs. Res. 

96, 1159, 1991. 

Discuss observed vs. adjusted for F'10.7 values, also 
URSI Series D 

For testing http: //nssdc.gsfc. nasa.gov/space/model/models/msis 

http://www.agu.org/joumaIs/ja/ja0212/2002.JA009430/ 
go to auxiliary material on the left side menu and open 
the tables-datasets.doc 

Other useful models http: / /nssdagsfe.nasa. gov/space / model / 

Exponential Atmosphere 
Solar Radiation Pressure 

5.2.5 Spacecraft Thrust 

5.3 Environment Modelling 

5.3.1 Celestial Body Ephemeris 

5.3.2 Analytic Ephemeris Model 

• For a new body, the user must input the central 
body by choosing from the 9 Planets or the sun. 

• the user must provide the epoch. 

• The user must provide the keplerian elements, in 
the central body centered, MJ2000Eq axis system. 


5.2.4 Solar Radiation Pressure 


a., 


,, CrA, 
-I'sr s 


TO, 


(5.71) 


where si is a unitized vector pointing from the spacecraft 
to the sun 

s = r s — r (-5.72) 

where r s is the Sun’s position vector and r is the space- 
crafts position vector. 
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Chapter 6 

Attitude 


The attitude of a spacecraft can be defined qualita- 
tively as how the spacecraft is oriented in inertial space, 
and how that orientation changes in time. GMAT has 
the ability to model the orientation and rate of rotation 
of a spacecraft using several different mathematical mod- 
els. Currently, GMAT assumes that a spacecraft is a rigid 
body. 

There are many ways to quantitatively describe the 
orientation and rate of rotation of a spacecraft, just like 
there are many ways we can quantitatively describe an 
orbit state. Let’s define any set of numbers that can 
uniquely define the spacecraft attitude as an attitude pa- 
rameterization. GMAT allows the users to use several 
common attitude parameterizations including quaternions, 
Euler angles, the Direction Cosine Matrix (DCM), Euler 
angle rates, and the angular velocity vector. Given an 
initial attitude state, GMAT can propagate the attitude 
using one of several kinematic attitude propagation mod- 
els. 

In this chapter, we discuss the attitude parameteriza- 
tions supported in GMAT, and how to convert, between 
the different types. We discuss the internal state param- 
eterization that GMAT uses. Next we investigate the 
types of attitude modes in GMAT and discuss in detail 
how GMAT propagates the spacecraft attitude in all of 
the Kinematic attitude modes. We conclude the chapter 
with a discussion of how GMAT converts between differ- 
ent attitude parameterizations. 


6.1 Attitude Propagation 

Given a set of initial conditions that define the attitude, 
GMAT can propagate the attitude using several methods. 
Currently, GMAT only supports kinematic attitude prop- 
agation. In Kinematic mode, the attitude is defined by 
describing the desired orientation with respect to other 
objects such as spacecraft or celestial bodies. With this 
information, GMAT can calculate the required attitude 
to satisfy the desired geometrical configuration. This sec- 
tion presents the different Kinematic attitude modes, and 
how GMAT calculates the attitude state in each mode. 


Let’s begin by looking at the internal attitude state rep- 
resentation and how the user can define initial conditions. 

6.1.1 Internal State Representation and 
Attitude Initial Conditions 

Certain attitude parameterizations are more useful for 
attitude propagation, while other attitude parameteriza- 
tions are more intuitive for providing attitude initial con- 
ditions or output. GMAT uses different internal param- 
eterizations of the attitude orientation depending upon 
the attitude mode. The type of parameterization is cho- 
sen to make the attitude propagation algorithms natural 
and convenient. For the kinematic modes, GMAT uses 
the DCM that represents the rotation from the inertial 
system to the body axes as the attitude orientation pa- 
rameterization. In the future, when 6 degree of freedom 
attitude modelling is implemented, GMAT will use the 
quaternion that represents the rotation from the inertial 
system to the body axes. GMAT uses the angular velocity 
of the body with respect to the inertial frame, expressed 
in the body frame, {^ib}bi as the rate portion of the 
state vector. 

For convenience, the user can choose a coordinate sys- 
tem in which to define the initial attitude state. Let’s 
call this system Ti. The user can define the initial at- 
titude orientation with respect to Ti using Euler angles, 
the DCM, or quaternions. The user can define the body 
rate with respect to T t by defining the angular velocity 
in Ti, {u ’ ib}x, or by defining the Euler angle rates. Note 
that not all attitude modes require these three pieces of 
information. The specific inputs for each attitude mode 
are discussed below, along with details about how atti- 
tude propagation is performed in each mode. 

6.1.2 Kinematic Attitude Propagation 

The Kinematic attitude mode allows a user to define a 
geometrical configuration based on the relative position 
of a spacecraft with respect to other spacecraft or celes- 
tial bodies, and with respect to different coordinate sys- 
tems. In Kinematic mode, GMAT does not integrate the 
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attitude equations of motion, but rather calculates the 
attitude based on the geometrical definition provided by 
the user. There are several Kinematic modes to choose 
from. The different modes allow the user to conveniently 
define the spacecraft attitude depending on the type of 
attitude profile needed for a specific mission. To begin, 
let’s look at how GMAT calculates the attitude state in 
the Coordinate System Fixed attitude mode (CSFixed). 

Coordinate System Fixed Mode 

In the CSFixed attitude mode, the user supplies two pieces 
of information. They first specify a coordinate system in 
which to fix the attitude, Ti. T, can be any of the de- 
fault coordinate systems or any user defined coordinate 
system. Secondly, the user specifies how the body axis 
system, Tb is oriented with respect to Ti by defining 
Ilfli or an equivalent parameterization. With this infor- 
mation, GMAT calculates the rotation from the inertial 
to the body axes and the angular velocity of the body 
with respect to the inertial frame, expressed in the body 
frame, {u IB } B . 

GMAT calculates the rotation matrix from Ti to Tb , 
Rsi, from the initial conditions provided by the user. For 
CSFixed mode, R B « is constant and is stored for use in 
the equations below. Knowing R B< , we can. calculate the 
rotation matrix from the inertial frame to the body frame, 
Rij/ , using the following equation 

R«/ — R.«iR>/ (6.1) 

R<f is the rotation matrix from Tj to Ti and GMAT 
knows how to calculate this matrix for all allowable Ti. 
For details on the calculation of this matrix for all coor- 
dinate systems in GMAT see Ch. 3. 

To calculate {cojb}b ■. we start from Eulers’ equation: 
R bi — — {w x /b}bF1.b/ (6.2) 

where 

/ 0 ~uj 3 v 2 \ 

{^ x ib}b — ^3 0 —oJi j (6.3) 

\—oJ2 uq 0 ) 

and {<jJib}b is the rotation of Tb with respect to Tj, 
expressed in Ti j. Solving Eq. 6.2 for we obtain 

{^ x ib}b -Rfl/Rgf (6.4) 

Taking the derivative of Eq. (6.1) with respect to time 
yields 

R bi = Rfl,Rff (6.5) 

because by definition, for the CSFixed mode, R B i = 0. 
Substituting Eq. (6.5) into Eq. (6.4) we obtain 


where Rg, is known from user input, and R bi is known 
from Eq. (6.1) . GMAT knows how to calculate R ( j for 
all allowable T% and details are contained in Ch. 3. 

In summary, in CSFixed mode, Eq.(6.1) is used to 
calculate R B /, and Eq. (6.6) is used to calculate 
If another attitude parameterization is required, GMAT 
uses the algorithms in Sec. 6.2 to transform from R B j 
and \ b to the required parameterization. Now let’s 
look at the spinning spacecraft mode. 

Spinning Spacecraft Mode 

In spinmiing spacecraft mode, GMAT propagates the at- 
titude by assuming the spin axis direction is fixed in in- 
ertial space. The spacecraft attitude at some time, t, is 
determined from the attitude initial conditions, the angu- 
lar velocity vector, and the elapsed time from the initial 
spacecraft epoch. Let’s take a closer look at the calcula- 
tions. 

In the spinning spacecraft mode, the user provides 
three pieces of information. They first choose a coordi- 
nate system, Ti, in which to define the initial conditions. 
Secondly, they define the initial orientation with respect 
to Ti by providing Rgj or an equivalent parameterzation 
that is then converted to the DCM. The user also pro- 
vides the angular velocity of the body axes with respect 
to the inertial axes expressed in T, {u>r B }i- 

To calculate R B /(f) where t. is an arbitrary epoch, we 
begin by calculating R b 0 i where R Bo j = R We 
calculate R b 0 i using 

Rfj 0 / = Rb?;R*j(G) (6.7) 

where R Bi comes from user provided data, and R;/(f„) 
is calculated by GMAT and is dependent upon Ti. See 
Ch. 3 for details on how GMAT calculates Ru for all 
allowable coordinate systems in GMAT. 

Before calculating Rg r ( t ) we must determine the spin 
axis in the body frame, {u>ib}b- The user provides {w/b }< 
In spinning mode we assume the spin axis direction is con- 
stant in inertial space and in the body frame so {w tb}b(£) 
= ( w /b)b(G) — {^/b}b- We can find the spin axis in 
the body frame using R B< as follows 

{w/bJb ~- RbjW/b } i (6.8) 

Once calculated, GMAT saves R, b„i and {u>/b}b for use 
in calculating the attitude orientation and rate at other 
epochs. 

GMAT calculates R B j(£) using the Euler axis/angle 
rotation algorithm in Sec. 6.2. The Euler axis is simply 
the unitized angular velocity vector or, 

{^ib}b 


{fy ,x ib}b — -Ri?iR,/R. B/ 


(6.6) 


a = 


^IB 


(6.9) 
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where 

UlB — ||{w/b}b|| 

(6.10) 

The Euler angle 6 

is calculated using 



rf>(t) — &iB(t — i c ) 

(6.11) 


6.2.2 Conversion: DCM to Quaternions 

Given: R 
Find: q, q 4 

Define following vector 


where t is the current epoch, and t„ is the spacecraft’s 
initial epoch. Let’s define the rotation matrix that results 
from the Euler axis/angle rotation using a and d>(t), as 
Rjjj3 (j (t). We can calculate R sr(f) using 

R-Br(^) = (6.12) 

To summarize, in spinning mode the user provides 
Rjji and GMAT assumes that that the spin 

axis direction is constant, and uses the Euler axis/angle 
method to propagate the attitude to find Re/. 

Now let’s look at how GMAT performs conversions 
between the different attitude parameterizations. 


6.2 Attitude Parameterizations 
and Conversions 


This section details how GMAT converts between differ- 
ent attitude parameterizations. For each conversion type, 
any singularities that may occur are addressed. The ori- 
entation parameterizations in GMAT include the DCM, 
Euler Angles, quaternions, and Euler axis/angle. The 
body rate parameterizations include Euler angle rates and 
angular velocity. We begin with the algorithm to trans- 
form from the quaternions to the DCM. 


v = [R n R -22 A 33 trace(R) ] (6.17) 

Define i m as the index of the maximum component of v 
if i m — 1 

l‘2vi m + 1 — trace) R)^ 

(6.18) 


if im — 2 


If hn — *1 


if L 


( 2v, 

m +1 

— trace) R)\ 


R \.2 

,/^21 



Rvz 



V 

A 23 

— A32 


( 

R 21 

+ A 12 

\ 

2v. 

m +1 

— trace(R) ] 


«J 3 

+ a 32 


V 

R .31 

- A] 3 


/ 

A 3 1 

+ A 13 

\ 


A3 2 

+ A23 


2v. 

m +1 

— trace(R) 

\ 

A], 2 

— A2] 

J 


( a 23 

-A 32 \ 



A31 

- An 



A12 

- A21 



\1 + traee(R)/ 


(6.19) 


( 6 . 20 ) 


( 6 . 21 ) 


We normalize q" using 


q = 


( 6 . 22 ) 


Finally, 

6.2.1 Conversion: Quaternions to DCM q =[q'r q' 2 q’z ] T (6-23) 

and 

Given: q, 94 g 4 ^ ^ (6.24) 

Find: R 


Name: ToCosineMatrix 


= (51 

92 9 a) 

T 

( 0 

-93 

92 

93 

0 

-91 

\-?2 

9 l 

0 


1 



q'i + q'i + <ii + q'i 


R = c [(ql - q T q)I.3 + 2qq T - 2q 4 q x ] 


6.2.3 Conversion: DCM to Euler 
Axis/ Angle 


(6-13) 

Given: R 



(6.14) 

Find: a, <p 



(6.15) 


/An A12 A 13 \ 
R == R n R 22 A 23 
\A 3 i A 3 2 R33J 

(6.25) 

(6-16) 


d> = cos -1 (trace(R) — l)'j 

(6.26) 



66 


CHAPTER : 6. ATTITUDE 


«23 - «S2\ 

ifei - Ri» 




If || sin ^|| < 10 , then we assume 


a = [ 1 0 0 ] T 

Note that if || sin <^>|| < 10 -14 then cos ^«1 and 
at a DCM of I3 . 


(6.27) 

(6.28) 
we arrive 


6.2.4 Conversion: Euler Axis/Angle to 

DCM 

Given: a, <b 
Find: R 


“«3 

0-2 \ 


0 

-oA 

(6.29) 

ttl 

° / 


<£)aa T 

— sin <p& x 

(6.30) 


= 0 c 3 s 3 0 1 0 


C! 0 


The approach is similar for the remaining 11 Euler 
angle sequences. Rather than derive the DCM matrices 
for the remaining 11 sequences, we present them in Table 
6 . 1 . 


6.2.6 Conversion: DCM to Euler Angles 

Given: Sequence order ( i.e. 123, 121, .... 313), R 
Find: 0\, 0 2 , 0 3 

Well give an example for a 321 rotation, and then 
present results for the remaining 11 Euler angle sequences. 
Examining, Eq. (6.35), we see that 

R 2 \ cos sin 0\ 

R,n cos 62 cos 0i 

From this we can see that 


(6.36) 


0 i = tan 


_i R21 
Ru 


(6.37) 


6.2.5 Conversion: Euler Angles to DCM 

Given: Sequence order ( i.e. 123, 121, .... 313), 0\, 62 , 63 
Find: R 

We’ll give an example for a 321 rotation, and then 
present results for the remaining 11 Euler angle sequences. 
First, let’s define R 3 (#i), R 2 (0 2 ), and Ri(0 3 ). 

( cos 0\ sin 0\ 0\ 

— sin cos(9i 0 I (6.31) 
0 0 l) 

( cos 02 0 —sin 02 \ 

0 10 (6.32) 

sin 02 0 cos 02 J 

A 0 0 \ 

Ri(0 s ) — I 0 cos 0.3 sin 0 3 I (6.33) 

\0 — sin 03 cos 03 / 

Now we can write 

R.321 ==:R 1 (0 3 )R 2 (02)R 3 (^) = 

(l 0 0 \ ( c 2 0 -s 2 \ / Cl A'l 0\ 


Further inspection of Eq. (6.35) shows us that 

0 2 = sin -1 R13 (6.38) 

At first glance, we may choose to calculate 0s using 0 3 = 
tan -1 (R23/R33). However, in the case that 0 2 — 90°, this 
would result in the indeterminate case, 03 — tan" 1 (R23/.R33) 
= tan 1 (0/0). An improved method, found in the AD EAS 
mathematical specifications document, is to determine 0 S 
using 

_i f ((31 sin0j 


0 3 — tan 


R21 sin 0 


— R32 COS 01 \ 

i + R-22 cos 01 ) 


(6.39) 


\0 -S3 \s 2 0 c 2 J \ 0 0 1 / 

(6.34) 

where ci — cos0i, .sq — sin0i etc. We can rewrite R.321 

/ C 2 Ci C 2 Sl — s 2 \ 

R32I = J —C3S1 + S3S0C. C3C1 + S 3 S 2 Sl S 3 C 2 I 
\ S3S1 + C3.S2C1 —S3 Cl -f c a s 2 S] C3C2/ 

(6.35) 


Substituting values from Eq. (6.35) into Eq. (6.39), and 
using abbreviated notation, we see that 

a . _i / si(s 3 si + C3S2C1) - ci (—s 3 ci 4 - C3S2S1) 

V Sl(C 3 Sl - S3S2C1) + Ci(C;)Ci + S 3 S 2 Sl) 

(6.40) 

Now, if 0 2 — 90°, and we substitute c 2 = 0 and s 2 — 1 
into the above equation, we see we get a determinate form. 
Results for all twelve Euler Sequences are shown in Table 
6.3. 

Note: All tan" 1 use a quadrant check ( equivalent to 
atan2 ) to make sure the the correct quadrant is chosen. 


6.2.7 Conversion: Angular Velocity to 
Euler Angles Rates 

Given: Sequence ( i.e. 123, 121, .... 313), 0 2 , 0 3 u> 
Find: 0], 0 2 , 0 3 


0 2 =S- 1 (0 2 ,0 3 )o; 

\0-J 


(6.41) 
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S -1 (0 2 , 0 3 ) is dependent upon the Euler sequence. Table 
6.2 contains the different expressions for S 1 {0 2 ,0a) for 
each of the 12 unique Euler sequences. 

Note: Each of the forms of S -! have a possible singu- 
larity due to the appearance of either sin 0 2 or cos 0 2 in 
the denominator. If GMAT encounters a singularity, an 
error message is thrown, and the zero vector is returned. 


6.2.8 Conversion: Euler Angles Rates to 
Angular Velocity 

Given: Sequence ( i.e. 123, 121, .... 313), 6 9 , 0 3 , (>\ , # 2 , 

k 

Find: u> 

(<>A 

w = S(0 2 ,0 3 ) (02 I (6-42) 

S(0 2 , 0 3 ) is dependent upon the Euler sequence. Table 6.2 
contains the different expressions for S _1 (02,0s) for each 
of the 12 unique Euler sequences. 


6.2.9 Conversion: Quaternions to Euler 
Angles 

Given: q, </ 4 , Euler Sequence 
Find: 0\ 1 O2, and 63 

There is not a direct transformation to convert from 
the quaternions to the Euler Angles, GMAT first converts 
from the quaternion to the DCM using the algorithm in 
Sec. 6.2.1. The DCM is then used to calculate the Eu- 
ler Angles for the given Euler angle sequence using the 
algorithm in Sec. 6.2.6. 

6.2.10 Conversion: Euler Angles to Quater- 
nions 

Given: fcfi, 0 2 - and 63, Euler Sequence 
Find: q, q A 

There is not a direct transformation to convert from 
Euler Angles to quaternions. GMAT first converts from 
the Euler Angles to the DCM using the algorithm in Sec. 
6.2.5. The DCM is then used to calculate the quaternions 
using the algorithm in Sec. 6.2.2. 


6.3 Appendix 1 
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Table 6.1: Rotation Matrices 
R'3(^3)R'2(^2)R-l(A) — 

R2(®3)Bfl(«a)R.i(«i) - 

R a (0. 3 )R 3 (^)R.2(^l) = 

R-3($3)Rl(# 2 )R2($l) = 

R2(0 3 )Ri(0 2 )R 3 (f/i) - 

Rl(^)R 2 (02)R3(^l) = 

R ] (0 3 ) R 2 (6> 3 )Ri (Ox ) = 

R i (# 3 ) R 3 (# 2 ) Ri (#i ) = 

R 2 (^)Ri(02)R2^l) = 

R2(^)R3(^2)R ? (^l) - 

R 3 (0 3 )Rl (02^3(00 = 


R 3 (^)R 2 (^)R3(^.) 
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Table 6.3: Computation of Euler Angles from DCM 


Euler Sequence 

Euler Angle Computations 



R 3 (0 3 )R2(02)Ri(0i) 

0\ = tan — ^ 32 /^ 33 ) 

0 2 = sin ” 1 (R 31 ) 

0 3 = tan" 

R.2(03)R3(02)Ri(0i) 

01 - tan” 1 (R 23 IR- 22 ) 

0 2 — sin” 1 (— fi 2 i) 

0 3 — tan" 

Rl (.03 )R 3 (0 2 )R 2 (01 ) 

0i =s tan“ 1 (— Ri$/Rii) 

0 2 — sin 1 (Ri 2 ) 

0 3 — tan" 

R 3 (0 3 )Ri(0 2 )R2(0i) 

0i — tan” 1 (fi 3 i/J?, 33 ) 

0 2 = sin _1 (— R 32 ) 

0 3 — tail' 

R. 2 (0 3 )Ri(0 2 )R3(0 1 ) 

0i — tan” 1 (— R 2 i/R 22 ) 

0 2 — sin -1 (R 2 s| 

03 — tail " 

R-i (0 3 ) R 2 (0 2 ) R 3 (0 i ) 

0i = tan” 1 (fi 12 /i?n) 

02 — sin _1 (— Ris) 

0 3 — tan' 

Rl(03)R 2 (02)Rl(0l) 

0i — tan'" 1 (fii 2 /(— firs)) 

0 2 — cos” 1 (fin) 

0 3 = tan’ 

Ri (0 3 )R 3 (0 2 )Ri (0i ) 

0i - tan- 1 (Ru/ (R- 12 )) 

09 — cos” 1 (fin) 

0 3 — tail' 

R-2 (03 ) Rl (02 ) R 2 (01 ) 

0i — t an" 1 ( /i 2 i / ( W 2 3 ))• 

0 2 — cos” 1 (fi 2 2 ) 

0 3 — tail 

R.2 (03 ) R'3 (02 ) R 2 (01 ) 

0i = tan 1 (fi 2 3/(— R 21 )) 

O 2 — cos ^(^ 22 ) 

0 3 — tan 

R 3 (03)Rl(0 2 )R3(0l) 

01 — tan” 1 (yf-3! /( R 32 }) 

O 2 — cos'” ^(^ 33 ) 

0 3 — tan 

R'3 (0 3 ) R 2 (0 2 ) R 3 (01 ) 

01 - tan” 1 ( fi .32 / (R 31 ) ) 

0 2 — COS ” 1 (fi 33 ) 

0 3 — tan 


( fii 3 sin 0 \ + R-12 cos 0 i \ 

\ Ha sin + Ha cos 9 y ) 

f Ri2sm0i — RiscosQi \ 

V —R32 sin 9 1 + Ha cos ~ 9 \ ) 

( fi 2 i sin <?i + fi 23 cos 0 L \ 

V Hix sin 9 t + A 33 cos 0 ] J 

( R -23 sin 9 \ — R21 cos 9 \ \ 
\—Ri .3 sin 6 1 + Rn cos 9 \ j 

( ^' A2 sin ^l + Hy, cosflx \ 

V R\2 sin + Rn cos 9 1 J 

/ fi 31 sin 9 \ — fi 32 cos 9 \ \ 
\ — fi 2 i sin 9 1 + fi 22 cos 9 \ ) 

( — -R33 sin 9 1 — R 3 2 cos 9 i \ 
\ fi 23 sin 9 \ + R 2 2 cos 9 i J 

( — Hq sin 9 1 + Rg cos 9 \ \ 
V — H\2 sin 9 1 + R .33 cos 9 \ ) 

f — R 33 sin 9 1 + R 3 i cos 9 \ \ 
\ — Ri'a sin + Rn cos0i J 

( — Ri 1 sin 0, — R \ 3 cos 0 i \ 
\ A. 31 sin 0i + JK33 cos 0 i J 

( H 2 2 sin 0i - R 2 i cos 0! \ 
\ l~t \2 sin 0 1 T R-x 1 cos 0 1 J 

f —R\ 1 sin 0i + Ryz cos0! \ 
\-R -21 sin 0i + R 22 cos 0i J 


Chapter 7 

Numerical Integrators 


7.1 Runge-Kutta Integrators 

The Runge-Kutta integration scheme is a single step method 
used to solve differential equations for n coupled variables 


estimate of the accuracy of the step; a second set of coef- 
ficients corresponding to this second integration scheme 
can be used to obtain a solution 


of the form 


+ 5t) = 


stages 

r [n) (t)+ 


c*k* {n) 


ar' 1 

dt 




(The superscript in this discussion refers to the vari- 
ables; hence f is the i th variable, and refers to all 
n variables.) The method takas an integration step, h, 
by breaking the interval into several stages (usually of 
smaller size) and calculating estimates of the integration 
result at each stage. The later stages use the results of 
the earlier stages. The cumulative effect of the integra- 
tion is an approximate total step St, accurate to a given 
order in the series expansion of the differential equation, 
for the state variables n(t + St) given the state n(t). 


3 = 1 


With care, the stage estimates kj and k* can be se- 
lected so that they are the same; in that case, the estimate 
of the error in the integration A- 71 ) can be written 


] stages 

A (n) = ! 


-c*) kj 


>» ! 


(The difference between the coefficients Cj — c* is the 
array of error estimate coefficients (ee) in this code.) 

Once the estimated error has been calculated, the size 


The time increment for a given stage is given as a 
multiple Oj of the total time step desired; thus for the 
i th stage the interval used for the calculation is ajSt: the 
estimate of the integrated state at this stage is given by 

i- 1 

= 5tf(t + aiSt, r (n) ( t ) + ^2 

3=1 


of the integration step can be adapted to a size more ap- 
propriate to the desired accuracy of the integration. If 
the step results in a solution that is not accurate enough, 
the step needs to be recalculated with a smaller step size. 
Labeling the desired accuracy a and the obtained accu- 
racy e (calculated, for instance, as the largest element 
of the array A), the new step used by the Runge-Kutta 
integrator is 


where b a contains a set of coefficients specific to the 
Runge-Kutta instance being calculated. Given the results 
of the stage calculations, the total integration step can 
be calculated using another set of coefficients, c.j and the 
formula 

stages 

r< n >(i + St) = r (n >(t) + J2 
3=1 

The error control for these propagators is implemented 
by comparing the results of two different orders of inte- 
gration. The difference between the t wo s teps provides an 


/ a \ i/(»» — i) 

Stnew — (7 St ■ j 

where m is the order of truncation of the series ex- 
pansion of the differential equations being solved. The 
factor a is a safety factor incorporated into the calcula- 
tion to avoid unnecessary iteration over attempted steps. 
Common practice is to set this factor to 0.9; that is the 
default value used in this implementation. 

Similarly, if the step taken does not result in the de- 
sired accuracy, you may want to increase the step size 
parameter for the next integration step. The new esti- 
mate for the desired stepsize is given by 
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51 


lew = vfik ( — ) 
\ € / 


l/(?n) 


set of first order differential equations. The algorithm is 

found at http: //chemical . caeds . eng.uml . edu/onlinec/white/m< 

or iri Bate, Mueller and White, pp. 415-417. 


Sometimes you do not want to increase the stepsize 
in this manner; for example, you may want to keep the 
maximum step taken at some fixed value. This implemen- 
tation provides a mechanism for specifying a maximum 
allowed step. 


The predictor step extrapolates the next state ?q + i of 
the variables using the the derivative information (/) at 
the current state and three previous states of the vari- 
ables, by applying the equation 


Sometimes it is convenient to request steps of a spec- 
ified size, regardless of the stepsize control algorithm or 
the calculation of the ’’best step” described above. This 
implementation accomplishes that task by taking multi- 
ple error controlled steps is necessary to step across the 
requested interval. 


r *j ^ J 4 - — 
'<+). '» • 24 


tofl-MfLi+Wfl-2-9fLs 


The corrector uses derivative information evaluated 
for this state, along with the derivative information at 
the original state and two preceding states, to tune this 
Both of these features are implemented using the booleanstate, giving the final, corrected state: 
flags described in the base class for the integrators. See 
the documentation for the Integrator (p. ??) class for 
more information about these flags. . /, 


't+i 


= r? + 24 [9/£ i + ~ 5/n-i + 1 51 - 


7.1.1 


Constructor & Destructor Documen- 
tation 


Bate, Mueller and White give the estimated accuracy 
of this solution to be 


RungeKutta.uRungeKutta (int st, int order) 


Provides the greatest relative error in the state vector. 

This method takes the state vector and calculates the 
error in each component. The error is then divided by 
the change in the component. The function returns the 
largest of the resulting relative errors. 

Override this method if you want a different error es- 
timate for the stepsize control. For example, we are using 


errori = 


( A j(t + dt) \ 

\r»(f + St) — i ~i(t) J 


Method used to fire the step refinement (the corrector 
phase). 

7.4 Bulirsch-Stoer 

7.5 Stopping Condition Algorithm 

7.6 Integrator Coefficients 


Another popular approach is to divide the estimated 
error A, by the norm of the corresponding 3-vector; for 
instance, divide the error in x by the magnitude of the 
displacement in position for the step. 


7.2 Prince-Dormand Integrators 

7.3 Adams Bashforth Moulton 


Implementation of the Adams-Bashford-Moulton Predictor- 
Corrector. 


This code implements a fourth-order Adams- Bash ford 
predictor / Adams-Moulton corrector pair to integrate a 
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Table 7.1: Prince- Dormand 45 Coefficients 
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Table 7.2: Prince-Dorm and 56 Coefficients (Warning: There is an error in the original source for these and we have 
not found the correct coefficients yeti!) 
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Table 7.3: Runge-Kutta-Fehlberg 56 Coefficients 


a, ; 




b{j 





0 

0 








3. 

6 

6 

0 







4 

4 

16 

0 






15 

75 

75 






2 

5 

S 

6 

0 





3 

6 

3 

2 





4 

8 

144 

-4 

16 

0 




5 

5 

26 

25 




i 

361 

18 

407 

11 

55 

0 



320 

5 

128 

80 

128 



0 

11 

640 

0 

11 

256 

11 

160 

11 

256 

0 

O' 


l 

93 

18 

803 

11 

99 

0 


0 

640 

5 

256 

160 

256 



7 

0 

1125 

9 

125 

0 

M 


C 3 

1408 

2816 

32 

768 

66 

66 

e j 

7 31 

1408 " 384 

0 

0 

0 

0 

5_ 

66 

5 

66 

5 

66 





Chapter 8 

Mathematics in GMAT Scripting 


8.1 

Basic Operators 

the same length, dotp is the scalar product. 

8.2 

Math Functions 


8.2.1 

max 

8.2.6 cross 

[maxX] 

= max(X) 

[crossp] = cross (vecl, vec2) 

X is an nxm array. maxX is a lxm row 7 vector contain- 
ing the maximum value in each column of X. 

The cross function calculates the cross product of 
two vectors, vecl and vec2 must both be vectors with 
the same length, crossp is the cross product. 


8.2.2 min 

[minX] = min(X) 

X is an nxm array. minX is a 1 xm array containing the 
minumuin value contained in each row of X. 

8.2.3 abs 

[absX] - abs(X) 

X is an nxm array. absX is a nxm array where each 
component is the absolute value of the corresponding 
component of X. 

8.2.4 mean 

[meanX] = mean(X) 

X is an nxm array. meanX is a lxm row vector con- 
taining the mean of each column of X. 


8.2.7 norm 

[normv] = norm(vec) 

The norm function calculates the 2-norm of a vector, 
vec must both be a vector, normv is the Toot-sum-square 
of the components of vec. 

8.2.8 det, 

[detX] = norm(X) 

The det function calculates the determinant of a ma- 
trix. X is an nxn array. detX is the determinant of X. 


8.2.5 dot 

[dotp] = dot (vecl, vec2) 

The dot function calculates the dot (scalar) product 
of two vectors, vecl. and vec2. must, both be vectors with 


8.2.9 inv 

[invX] = inv(X) 

The inv function returns the inverse of a matrix. X 
must be a square matrix. 
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8.2.10 eig 

8.2.11 sin, cos, tan 

8.2.12 asin, acos, atan, atan2 

8.2.13 sinh, cosh, tanh 

8.2.14 asinh, acosh, atanh 

8.2.15 transpose 

8.2.16 DegToRad 

8.2.17 RadToDeg 

8.2.18 log 

8.2.19 loglO 

8.2.20 exp 


8.2.21 


sqrt 




Chapter 9 

Solvers 

9.1 Differential Correction 

9.2 Broyden’s Method 

9.3 Newton’s Method 

9.4 Matlab fmincon 


The user first creates a solver and names it. An example 
is 


Create fminconGpt iniizer SPQfmincon 

The user creates an optimization sequence by issuing 
an optimize command, followed by the name of the opti- 
mizer to use 

Optimize SQPfmincon 

EndOptimize 

9.5 The Vary Command 

The user defines the independent variables by the vary 
command, 


Table 9.1: Available Commands in an fmincon Loop 


Value 

Command 

Xt 

Vary 

Upper Bound on X, 

Vary 

Lower Bound on X* 

Vary 

Nondimensionalization 

Vary 

Factor 1 

Nondimensionalization 

Vary 

Factor 2 

Nonlinear constraint 

NonLinearConstraimt 

function 

Linear constraint bin c- 

Linear Constraint 

tion 

Cost Function 

OptimizerName .Cost = 
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Chapter 10 

Event Functions 


Event functions in GMAT allow a user to determine 
when different types of events occur such as station con- 
tacts, eclipse events, or spacecraft-to-spacecraft line of 
sight. In general, these events are dependent upon the 
orbit state and other time dependent parameters, and 
therefore can only be determined during or after orbit 
propagation. The implementation of Event Functions re- 
quires GMAT to find the roots of a parametric function of 
time. The roots of the parametric equation are the event 
times. 

In this chapter, we’ll look at how GMAT calculates 
the roots of Event Functions. For the implementation 
in GMAT, this includes two subproblems. The first is 
determining if a root, has occurred during a propagation 
step. The second, is determining the numerical value of 
the root. There is a trade between robustness and perfor- 
mance, so well look at several options provided in GMAT 
that allow the user to select between a more robust yet 
slower approach, or a fast but less robust approach. Let’s 
begin my looking at the mathematical definition of an 
event function in GMAT. 

10.1 Event Function Mathemati- 
cal Definition 

An Event Function in GMAT has three outputs. The 
general form of an Event Function is 

[ f»d,p ] = F(f,x(t), C) (10.1) 

where t is the current time, x(f) is a vector of time de- 
pendent parameters such as spacecraft states, and C is a 
vector of constants, f is vector of function values at t, d 
is a vector describing the or sign change we wish to track 
that occurs at the root, and p is a vector that tells GMAT 
whether a root is possible or not. Let’s talk about some 
of the output variables in more detail. 

For efficiency and convenience, the user can calculate 
several different function values, /, inside of a single event 
function, F. This is useful when several functions require 
similar yet expensive calculations. GMAT allows the user 


to pass back a vector of function values in the; output 
parameter f, where the components of f are simply the 
values of the different functions /, or 

f = [ /i(f.x(t), C), f 2 (t,x(t), C C) } T 

(10-2) 

The output parameter d allows the user to define which 
type of roots for GMAT to calculate. For example, in 
some cases we might only be interested in roots that oc- 
cur when the function changes from a negative value to a 
positive value. In other situations we may only be inter- 
ested in roots that occur when the fimction passes from 
positive to negative. Finally, we may be interested in both 
types of roots, d is a vector that has the same number of 
elements as f, and the first element of d corresponds to 
the first element of f and so on. Table 10.1 summarizes 
the allowable choices for components of d and the action 
GMAT will take depending upon the selection. 


Table 10.1; Allowable Values for d in Event Function 
Output 


Value 

Action 

d = 1 

Find roots when the fimction is moving in 
the positive direction. 

d = -l 

Find roots when the fimction is moving in 
the negative direction. 

a. 

II 

to 

Find both types of roots 


The last output variable in the Event Fimction output, 
p, is a flag that allows the user to tell. GMAT whether or 
not a root is possible. If a component of p is zero, then 
GMAT will not attempt to try to find a root of the corre- 
sponding function in f. This flag is included to improve 
the efficiency of the algorithm. It is often possible to per- 
form a few simple calculations to determine if a root is 
possible or not. For example, let’s assume an. event func- 
tion is written to track Earth shadow crossings and that 
the fimction is positive when a spacecraft is not in Earth’s 
shadow, and negative when it is in Earth shadow. It is a 
relatively simple calculation to determine if a spacecraft 
is on the day-side of Earth by taking the dot product of 
the Sun vector and the spacecraft’s position vector. If the 
quantity is positive, there is no need to continue calculat- 
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ing the actual function value. 

Now that we’ve looked at the definitions of the inputs 
and outputs of an Event Function, let’s look at some dif- 
ferent approaches to finding the roots of an event func- 
tion. 

10.2 Issues in Locating Zero Cross- 
ings 

Before discussing the practical issues in finding roots of 
Event Functions, let's take a look at a hypothetical func- 
tion to illustrate some of the issues that must be ad- 
dressed. Figure 10.1 shows a sample event function. The 
smooth line represents the locus of points of the function 
itself, and the large “X” marks represent the function 
values at the integration time steps. The smaller tick 
marks indicate the function values at the internal inte- 
grator stages, which may be available if we use a dense 
output numerical integrator. 

In general, we don’t have continuous time expressions 
for the inputs to event functions. We only know the in- 
puts to Event Functions at discrete points in time, so we 
only know the Event Function values at discrete points 
in time. Since these discrete times come from the nu- 
merical integration of a differential equation, we can only 
calculate Event Function values at the integration time 
steps, or at the internal stages if the information is avail- 
able. This fact can cause a significant problem because an 
Event Function may vary rapidly and the discrete times 
at which we know the Event Function may not give an 
accurate picture of the function. 

Let’s consider a few ways in which we can determine 
if a root has occurred, given a set of times and Event 
Function values. The most obvious method is to sim- 
ply look for sign changes in the function values. If the 
function changes sign, then we know we have bracketed 
a root. This approach will incorrectly conclude that a 
zero crossing did not occur, if there is an even number 
of zero crossings between two function values. Another 
approach is to fit a polynomial through the data, and see 
if the polynomial has any real roots. While this approach 
may be more accurate than looking for sign changes in 
some cases, it still does not guarantee a zero crossing is 
missed. A third approach might be to force the integra- 
tor to select step sizes based on the rate of change of 
the Event Function. We will not investigate t his method 
further here though. 

In short, there is no way to guarantee that a root 
crossing is missed. However, by having an understanding 
of the Event Function, having control over the maximum 
integration step size, and having access to the internal 
integrator stages when using a dense output integrator. 


we can do a acceptably good job determining when zero- 
crossings occur. Once a zero-crossing is identified, there 
are well known ways to calculate the actual root value. 

F — F(f,x(£),C) (10.3) 

We want to find all t such that 

F(f, x(f), C) = 0 (10.4) 


10.3 Root Finding Options in GMAT 

In implementing a root finding approach, we need to bal- 
ance accuracy and the need to find every root, with speed 
and performance. One way to do this is to allow the user 
to select between different approaches depending upon 
the accuracy needed for a particular application. The user 
has several controls to tell GMAT how to determine if a 
zero crossing has occurred, and how to calculate the nu- 
merical value of a root if one has been detected. Let’s look 
at the choices implemented in GMAT, and discuss some 
options that can be included if a more robust method is 
required. 

The first group of controls available to the user are 
related to how or if GMAT tries to determine if a root 
has occurred. The user can provide a flag in the output 
of an Event Function that tells GMAT whether it is pos- 
sible that a root has occurred during the last integration 
step. Tliis flag is notated as p and Is discussed in section 
10.1. If an element of p is zero, then GMAT will not use 
more sophisticated and therefore more computationally 
intensive methods to determine if a zero crossing for the 
particular component of the event function has occurred, 
p can be either zero or one, and can change value during 
propagation. If p changes from zero to one, GMAT be- 
gins using a root checking method specified by the user 
to determine if a zero crossi ng has occurred, and begins 
storing function data in case it is needed to interpolate a 
root location. 

The second control that determines if a zero crossing 
has occurred is called RootCheekMethod in the GMAT 
script language. There are several RootCheekMethod 
options available and the user can currently select be- 
tween FunctionSignChange and PolynomialFit. If the 
user selects FunctionSignChange, then GMAT looks for 
sign changes in the function output to determine if a zero 
crossing has occurred. If the user selects PolynomialFit, 
then GMAT fits a polynomial to the Event Function data, 
and checks to see if the polynomial has any real roots. 

If the polynomial has real roots, then a zero crossing 
has occurred. The type of polynomial GMAT uses in 
RootCheekMethod is the same as it uses in RootFindingMethod 
and is discussed in more detail below. 
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Figure 10.1: Sample Event Function Output 


If a zero crossing is detected, there are many ways 
to determine the numerical value of the root. The user 
can select between the different methods by using the 
RootSolvingMethod option. The two methods currently 
implemented in GMAT are called Quadrat icPolynomial 
and CubicSpline in the GMAT script language. As the 
name suggests, if the user selects QuadraticPolynomial, 
then GMAT uses the last three function values to cre- 
ate a quadratic polynomial. Then, the quadratic equa- 
tion is used to determine the root locations. Similarly, if 
the user selects CubicSpline, GMAT constructs a cubic 
spline and then uses interpolation to find the root value. 

Allowing the options above requires that care is taken 
in designing an algorithm to track events. Tn the next 
section we discuss some of the issues that must be ad- 
dressed in the Event Function algorithm, and present a 
flow chart that, describes the algorithm in detail 


10.4 Algorithm for Event Functions 


Table 10.2: Variables in Event Function Algorithm 
Variable Definition 

Number of data points required to use the 
requested RootSolvingMethod option 
n c Number of data points required to use the 

requested RootCheckMethod option 
f A vector of function values provided by 

the user defined Event Function 
N The length of f , which is the number func- 

tion values contained in the output of a 
user defined Event Function, 
d A vector of flags (length N) that defines 

which type of roots to track, (negative to 
positive, positive to negative, or both) 
p A vector of flags (length N ) that tells 

whether or not a zero crossing is possi- 
ble. A component of p is one if a root is 
possible, otherwise it is zero. 

Startup A vector of flags of lenght N. The compo- 
nents of Startup correspond to the com- 
ponents of f. A component of Startup is 
one, if there is less than max(n r ,n c ) data 
points saved for use in root finding. Oth- 
erwise, a component, of Startup is zero. 
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Loop over the number of elements in Event Function output 
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10.5 


Appendix 1: Root Finding Al- 
gorithms 


where the values for cc, are known from the inputs. 

To calculate the coefficients a*, b t , a, and a!* we start 
by calculating the eight quantities 


10.6 Quadratic Polynomial 


We are given three data points defined by a vector of 
independent variables 


hi 

Ai 


x i -f 1 x i 

V i \ 1 - 


(10.15) 

(10.16) 


Next we solve the following system of linear equations 


x = [ xi X'i xs f r (10.5) 

and a vector of corresponding dependent variables 

y :=: [ Vi 2/2 2/3 ] T (10.6) 


we wish to find a quadratic polynomial that fits the data 
such that 

y = Ax 2 + Bx + C (10.7) 

We begin by forming the system of linear equations 


( x\ Xi \ 
x 2 x 2 1 
V 4 x 2 1 


A 

B 

C 



We can solve for the coefficients using 


2/i 

Xa 
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2/2 
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x 2 

^2 
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T 2 



X \ 

Vi 

i. 
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2/2 

•i 

„.2 


q 

a: 3 

Vs 

1 

x{ 

Xi 

1 

r 2 

x 2 

X2 

1 

T 2 

x 3 

X3 

1 


C = y-< — Ax 2 — Bx 3 


( 10 . 8 ) 


(10.9) 


( 10 . 10 ) 


( 10 . 11 ) 


10.7 Cubic Spline (Not-a-Knot) 


We are given five data points defined by a vector of inde- 
pendent variables 


AS = B (10.17) 


where the components of 

An 

— 2/7,3 + hi 

Av2 

— 2 /i-j -f - 112 

'4 13 

= 0 

A21 

= 0 

cm 

CM 

= h- 3 -f 2 h A 

d 2 3 

— 2/7,3 H- /14 

A31 

hi 

hi + h 2 

A32 

hih 2 
(hi + h 2 ) 

A 33 

h l 

/13 -j- /i.4 

the components of B axe 


/in :=:: 

#21 = 
B 3 i = 


are given by 

(10-18) 

(10.19) 

( 10 . 20 ) 
( 10 . 21 ) 
( 10 . 22 ) 

(10.23) 

(10.24) 

2 (h 2 + fca) + ,, h:>h4 . (1 0.25) 
(h 3 + h A ) 

(10.26) 


6 (A 2 - Ai) (10.27) 

6 (A 4 - A s ) (10.28) 

6 (A 3 -A 2 ) (10.29) 


and S — [ Si S 3 S 3 ] T . We can solve for the components 
of S using Cramers Rule as follows 


Bn 

B21 

Bsi 

A 12 
A 22 
A32 

A 13 

A23 

A33 


|A| 


An 

Hu 

Ai;i 

-A -21 

i£>i 

^23 

Asi 

B:n 

/I33 

1 AI 


(10.30) 


(10.31) 


X — [ Xi X2 X3 x<< \ l (10.12) 

and a vector of corresponding dependent variables 

y = [ 2/1 2/2 2/3 2/4 2/5 ] T (10.13) 

we wish to find the four cubic polynomials, i = 1,2, 3,4, 
such that 


a _ B 31 — A 31 S 1 — A 32 S 3 
— • 

A 33 

Now we can calculate S 2 and S A using 


S 2 — 


h 2 S 1 + ES, 
h\ + h 2 


h A S 3 + h 3 S- A 

h-3 + h A 


(10.32) 


(10.33) 


(10.34) 


Pi — mix — I,) 3 + bj(x - Xi) 2 + Ci(x — Xi) + di (10.14) 
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Finally, the coefficients for the i tri cubic polynomial 
are given by 



c. , c. 

6 hi 

Si 

2 

f/i+l Vi 2 hi Si 4" hiSi-±-\ 

hi 6 

Vi 


( 10 . 35 ) 

( 10 . 36 ) 

( 10 . 37 ) 

( 10 . 38 ) 


o 
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Chapter 11 

Spacecraft Model 


11.1 

Orbit 

11.2 

Attitude 

11.3 

Ballistics and Mass 

11.4 

Actuators 

11.4.1 

Thrust and Impulse Models 


GMAT uses polynomial expressions for the thrust and 
specific impulse imparted to the spacecraft by thrusters 
attached to the spacecraft. Both thrust and specific im- 
pulse are expressed as functions of pressure and temper- 
ature. The pressure and temperature are values obtained 
from fuel tanks containing the fuel. All measurements 
in GMAT are expressed in metric units. The thrust, in 
Newtons, applied by a spacecraft engine is given by 


P T { P,T) = {G 1 + C 2 P + C 3 P 2 + C 4 P c * + C 6 P c 7 + C 8 P Cs } +C 10 Cfj 12P } 


l+Ci X +CuP 


Pressures are expressed in kilopascals, and tempera- 
tures in degrees centigrade. The coefficients Cl - C14 
are set by the user. Each coefficient is expressed in units 
commiserate with the final expression in Newtons; for ex- 
ample, Cl is expressed in Newtons, C2 in Newtons per 
kilopascal, and so forth. 

Specific Impulse, measured in m/s (or, equivalently, 
Newton Seconds/kilogram) is expressed using a similar 
equation: 


1* P (P,T) - {AT + K 2 P + K s P 2 + K 4 P Kr ’ + K g P k ~ + K 8 P Ka + K 10 K*™ P } 


T 


T, 


ref 


H -/Cm t-lfiuP 
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Table 11.1 shows the default values for these coeffi- 
cients. 


11.5 Sensors 

11.6 Tanks 


Mass is depleted from the fuel tanks using equation 


dm Fr 
dt I $p(J 


( 11 . 1 ) 


where the thrust and specific impulse are given by ?? and 
11.4.1. 


This mass depletion is integrated along with the other 
parameters during propagation. 

The tank model in GMAT manages the fuel mass and 
the input variables for the thrust and specific impulse 
polynomials. The tank can be run in either a blow-down 
or pressure regulated mode. In pressure regulated mode, 
the pressure in the polynomial is held at a fixed value. In 
blow-down mode, the pressure decreases as fuel is used, 
following the ideal gas law: 


P V = nRT (11.2) 

In GMAT’s blow-down model, the temperature T and 
the number of pressurant molecules n in the tank are held 
constant, so the right side of this equation is constant. 
The gas volume available in the tank grows as fuel is con- 
sumed, and the pressure decreases accordingly. The gas 
volume V G in the tank is computed from the total tank 
volume, Vt, the mass of the fuel, Mp, and the density of 
the fuel, p: 


Table 11.2: Default Riel Tank Parameters 


Parameter 

Default Value 

Units 

FuelMass 

756 

kg 

Pressure 

1500 

kPa 

Temperature 

20 

C 

Reffemperature 

20 

C 

Volume 

0.75 

m 3 

RielDensity 

1260 

kg/m 3 

PressureRegulated 

true 



V a = V T -^ (11.3) 

Table 2 shown the default values for the tank param- 
eters. 



11.6. TANKS 


Table 11.1: Default Thrust and Specific Impulse Coefficients 


Thrust Coefficient 

Default 

Units 

Impulse Coefficient 

Default 

Units 

Ci 

500 

N 

A'i 

2150 

m/s 

Cl 

0 

N/kPa 

Kz 

0 

m/ (s • kPa) 

Ci 

0 

N/kPa 2 

Ki 

0 

m/(s • kPa 2 ) 

C4 

0 

N/kPa c ® 

a 4 

0 

m/ (s ■ kPa Kr ') 

c 5 

0 

none 

a 5 

0 

none 

Ck 

0 

N/kPa c? 

Ac 

0 

m/(s • kPa K7 ) 

C 7 

0 

none 

A 7 

0 

none 

c s 

0 

N/kPa c » 

As 

0 

m/(s • kPa K9 ) 

Q> 

0 

none 

a 9 

0 

none 

C w 

0 

N 

A10 

0 

m/ s 

C n 

1 

none 

A n 

1 

none 

c 12 

0 

l/'kPa 

A12 

0 

1/kPa. 

Cis 

0 

none 

a 13 

0 

none 

C u 

0 

1/kPa 

A m 

0 

1/kPa 
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Chapter 12 

MathSpecAppendices 


Vector Identities 


5a T a n T da. 
dx ‘ a dx 

(12.1) 

da d , T > 1/2 a. 1 da. 

dx dx ^ a dx 

(12.2) 

da ~ 1 d , T > -i/2 a r 5a 

<9a; ' a a a 3 dx 

(12.3) 

dr a ^ 1 5a aa T 5a 

dx \a 3 J o 3 dx " a 5 5a; 

(12.4) 
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CHAPTER 3. COORDINATE SYSTEMS 


• y- axis: Completes the right-handed set, 

• z-axis: Normal to the equatorial plane. 


Finally, for the Earth, the Equator axis system a true 
of date equator system and is calculated using the algo- 
rithm described in Sec. 3.4.5. 


The Equator system in GMAT is a true equator of date 
axis system. The equatorial coordinate system is defined 
only for celestial bodies. For a particular body, the equa- 
torial system is defined by the bodies equatorial plane 
and its intersection with the ecliptic plane, at the current 
epoch. The Earth and Moon have highly accurate mod- 
els for their equatorial systems and and are treated at the 
end of this section. For the remaining bodies in the solar 
system, the equatorial coordinate system is calculated in 
GMAT using data published bjr the International Astro- 
nomical Union (IAU). 4 The IAU publishes data that gives 
the spin axis direction and prime meridian location of all 
the planets and major moons as a function of time. For 
the Earth, GMAT uses FKo reduction for the Equator 
system. For the Moon, GMAT can use either the TAU 
data, or Eider angles provided in the JPL DE405 files. 

Let’s look more closely at the data provided by the 
IAU. Figure 3.5 contains an illustration of the three vari- 
ables, a G , S a , and W , that are used to define a body’s spin 
axis and prime meridian location w/r/t MJ2000Eq. a B 
and 5 0 are used to define a body’s spin axis direction. W 
is the body’s sidereal time. The equations for a 0 , S a , and 
W for the nine planets and the Earth’s moon are found 
in Tables 3.1 and 3.2. From inspection of Fig. 3.5 we see 
that 


R j 2k ,i = R[ (90° + a 0 )Rf (90° - S a ) (3.52) 

a 0 and 8 0 vary slowly with time, so we can assume the 
derivative of R« for the Equator system is the zero ma- 
trix. 

/ 0,0 0.0 0 . 0 \ 

R, hk ,i = 0.0 0.0 0.0 (3.53) 

\ 0.0 0.0 0 . 0 j 

If the user chooses to use the D E405 files to determine 
the Moon’s orientation, then GMAT gets a set of Euler 
angles and rates from the DE-405 files. We then use the 
following equations to determine and R./-,..i- 

= R3(0i) T Ri(02) T (3.54) 

- R 3 (0i) t R.i (0 2 ) + R 3 (0i)Rf (02)(3.55) 


3.4.2 MJ2000 Ecliptic (MJ2000Ec) 

The MJ2000 Ecliptic axis system is defined as follows: 

• a:- axis: Along the line formed by the intersection 
of the Earth’s mean equator and the mean ecliptic 
plane, at the J2000 epoch. The axis points in the 
direction of the first point of Aries. 

• y- axis: Completes the right-handed set. 

• z-axis; Normal to the Earth’s mean equatorial plane 
at the J2000 Epoch. 


The matrix to rotate from MJ2000 Ecliptic (MJ2000Ec) 
to MJ2000 Equatorial (MJ2000Eq) is a rotation about the 
x-axis through the obliquity of the ecliptic at the J2000 
epoch which is 23.439291°: 


/l.O 0.0 0.0 \ 

R - 0.0 0.91748206 -0.397777156 

\0.0 0.39777716 0.9174820621 


(3.58) 


GMAT uses more significant digits than included here. 
The rotation matrix is constant by definition so its time 
derivative is identically the zero matrix. 


/ 0.0 0.0 0 . 0 \ 

R= 0.0 0.0 0.0 (3.59) 

\ 0.0 0.0 0 . 0 / 


3.4.3 True of Epoch Equator (TOEEq) 


The True of Epoch Equator axis system is defined as fol- 
lows: 

• x-axis: Along the true equinox at the chosen epoch. 
The axis points in the direction of the first point of 
Aries. 

• j'-axis: Completes the right-handed set. 


where 




and 


R 3 (0i) 


0.0 

0.0 

0.0 \ 

0.0 

—02 sin 02 

02 COS 02 

0.0 

—62 COS 02 

—02 sin 02 J 


— 0 i sin 

0 i 

0\ COS 0\ 

0.0 

—01 cos 

0 i 

— 0 \ sin 9\ 

0.0 

0.0 


0.0 

0.0 


• z-axis: Normal to the Earth’s true equatorial plane 
at the chosen Epoch. 

(3.56) The TOEEq axis system is an intermediate system in 
FK5 reduction. Rj; and Rj.; for the TOEEq system are 
calculated using the following equations 

R/, - N T (t c )P T {t 0 ) (3.60) 

(3.57) where t„ is the epoch defined in the coordinate system de- 
scription provided by the user in the epoch field. Hence, 



